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SUMMARY 


During the past several months and since the inception of the Flight 
Data Statistics Office, a growing trend to analyze and process Earth 
Resources Flight Data has prevailed. In view of the growing interest 
at Marshall Space Flight Center to organize and establish capabilities 
to perform research studies in this area, endeavors in the areas of 
interpretation, analysis, and development of algorithms have provided 
the necessary computational programming tools for data processing and 
data handling and analysis. Algorithms that have been developed thus 
far, are adequate and have been proven successful for several pre- 
liminary and fundamental applications such as software interfacing capa- 
bilities, probability distributions, grey level print plotting, contour 
plotting, isometric data displays, joint probability distributions, 
boundary mapping, channel registration and ground scene classification. 
This report is written in two sections. Section I consists of the 
algorithms that have been developed individually under the existing 
contract and section II is a description of an Earth Resources Flight 
r, (ERFDP) , which handles and processes earth resources 
users control. 
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FOREWORD 


This final technical report, entitled "Correlation 
Signatures of Wet Soils and Snows", presents the results of 
18 months programming and algorithm development under the 
auspices of the National Aeronautics and Space Administration 
(NASA), with Mr. Robert Jayroe acting COR. This work was 
performed under contract NAS8-26797 (IITRI Project J6243) . 
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CORRELATION SIGNATURES OF WET SOILS AND SNOWS 

M. R. Phillips 
IIT Research Institute 

SUMMARY 

During the past several months and since the inception of 
the Flight Data Statistics Office, a growing trend to analyze 
and process Earth Resources Flight Data has prevailed. In view 
of the growing interest at Marshall Space Flight Center to or- 
ganize and establish capabilities to perform research studies 
in this area, endeavors in the areas of interpretation, analysis, 
and development of algorithms have provided the necessary compu- 
tational programming tools for data processing and data handling 
and analysis. Algorithms that have been developed thus far, are 
adequate and have been proven successful for several preliminary 
and fundamental applications such as software interfacing capa- 
bilities, probability distributions, grey level print plotting, 
contour plotting, isometric data displays, joint probability 
distributions, boundary mapping, channel registration and ground 
scene classification. This report is written in two sections. 
Section I consists of the algorithms that have been developed 
individually under the existing contract and section II is a 
description of an Earth Resources Flight Data Processor, (ERFDP) , 
which handles and processes earth resources data under a users 
control. 

SECTION I INTRODUCTION 

In this section a description of each algorithm developed 
under the existing contract will be presented. The algorithms 
were developed as building blocks toward an automatic data pro- 
cessor for processing earth resources flight data. Each algorithm 
has been demonstrated to be compatible in an overall automatic 
processing environment. These algorithms were developed using 
an IBM 7094 computer with 32K available core storage. 
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1. PROBABILITY DISTRIBUTION 

A probability density function was programmed, and the earth 
resources data processor which is available upon request as an 
option selection is included. This option gives a probability 
density function for any number of channels, not to exceed 12, 
of the multispec tral data. This option is used when calculating 
the slicing intervals of the probability distribution to provide 
grey level mapping of any ground scene image. An estimate of 
the probability density function is obtained digitally by divid- 
ing the range of the data into a desired number of class inter- 
vals d^. The probability of occurrence of a value in the 
class interval is then given by 

Pi<X) « N..(X)/N, 

where N^(X) is the number of values that occur within the range 
of d i <( X d^ + i and 

k-1 

n =Z n i 

i=0 

is the total size of the population where K = number of class 
intervals . 


A. Computer Program 

The computer program calculates the probability density of 
any channel of a multispec tral scanner consisting of 12 channels 
of data. The distribution table is then printed out for all 
channels requested, over all the data samples (see table I). 

The program also passes to the grey level module the probability 
distribution of each channel for grey level mapping. 
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TABLE I 

PROBABILITY DISTRIBUTION 


AMPLITUDE 

-1.5 

-1.4 

-1.3 

- 1.2 

- 1.1 

- 1.0 

- .9 

- .8 


1.2 

1.3 

1.4 

1.5 


CH 1 

CH 2 

CH 3 

CH 4 

CH 5 

CH 6 

CH 7 

CH 8 

CH 9 

CH 10 

CH 11 

CH 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2 

0 

2 

0 

0 

0 

0 

0 

1 

0 

0 

0 

3 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

9 

0 

5 

0 

0 

1 

0 

0 

1 

0 

0 

0 

21 

0 

7 

0 

0 . 

0 

1 

1 

8 

0 

5 

0 

29 

0 

15 

1 

0 

1 

0 

0 

22 

0 

4 

0 

35 

0 

17 

1 

0 

3 

2 

4 

31 

0 

8 

1 

39 

0 

34 

2 

3 

9 

7 

8 

47 

0 

12 

0 

48 

0 

39 

8 . 

8 

18 

20 

16 

60 

0 

24 

7 

55 

1 

47 

12 

12 

22 

25 

30 

55 

0 

21 

11 

66 

0 

51 

22 

22 

38 

40 

40 

40 

0 

38 

5 

58 

2 

55 

31 

41 

47 

48 

52 

27 

0 

45 

12 

49 

12 

49 

49 

70 

50 

55 

60 

11 

0 

57 

22 

42 

18 

44 

66 

82 

54 

62 

72 

6 

5 

59 

37 

33 

20 

39 

72 

96 

47 

51 

55 

5 

1 

62 

44 

25 

35 

33 

62 

79 

41 

45 

52 

0 

8 

58 

48 

19 

30 

25 

40 

62 

35 

39 

31 

1 

16 

51 

49 

12 

41 

18 

22 

38 

25 

25 

20 

0 

27 

47 

56 

6 

54 

7 

12 

12 

17 

16 

15 

0 

32 

35 

52 

1 

61 

8 

8 

10 

12 

8 

7 

0 

39 

25 

40 

0 

55 

2 

2 

7 

8 

5 

3 

0 

47 

22 

32 

1 

45 

1 

0 

2 

0 

1 

5 

0 

52 

8 

24 

0 

40 

0 

1 

1 

2 

0 

1 

0 

65 

5 

16 

0 

35 

1 

0 

0 

0 

0 

0 

0 

72 

7 

8 

0 

28 

0 

0 

0 

0 

0 

0 

0 

58 

3 

7 

0 

20 

0 

0 

0 

0 

0 

0 

0 

46 

1 

1 

0 

12 

0 

0 

0 

0 

0 

0 

0 

39 

0 

0 

0 

7 

0 

0 

0 

0 

0 

0 

0 

22 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

7 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


NUMBER OF OCCURRENCES 



NCH 

NSPS 

NSCANS 

NSTART 

NSTOP 

NBTLG 

MODE 

ITYPE 

MSFC 

LTN 

NSKIP 

NCRE 

XMAX 

XMIN 

NOCHS 

NWHICH 


B. Data Problem Parameters 

Number of channels or spectral bands on input tape 

Number of samples per record or resolution ele- 
ments across one scan 

Number of scans to process 

Starting resolution element 

Stopping resolution element 

Bit length of input data word 

Signifying FORTRAN or non- FORTRAN input tape 

Type of input data; Floating point or fixed point 

MSFC scanner format option 

Logical unit to load input tape 

Number of initial data records to skip before 
processing 

Data incrementation 

Maximum value in data set 

Minimum value in data set 

Number of channels to calculate probability density 
function 

Channel selection for probability distribution 
calculation 


2. GREY LEVEL MAPPING 

In order to preview earth resources data and obtain quick- 
look information, a grey level mapping program was written to 
in the ERFDP . This gives a pictorial display of the 
ground scene image quantized to 10 different levels. Characters 
are selected to represent varying shades of grey. The levels 
are calculated by slicing the probability distribution table of 
selected channels into 10 separate cells (see Figure 1). The 
number of occurrences in each cell is equally distributed over 
the probability distribution. This is done by calculating 
equal areas for each cell using the trapezoidal rule method. 

Each of the 10 areas would represent one cell of the probability 
distribution. Each resolution element is compared with all the 
cells, and the cell in which it falls is assigned the respective 
alphanumeric character. 
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ability Distribution Curve 


A. Computer Program 

The computer program was written as a module to be included 
in the ERFDP. The logic provides communication with the prob- 
ability distribution program, since information calculated in 
that module is used in grey level mapping. The data are input 
to the program module and each resolution element input to the 
module is compared with each cell unitl it falls within the 
quantized level of that cell. Each resolution element, x^, is 
examined such that 


SLICE (j-1) < x ± SLICE (j) 

where j = 1,2,3 ..., 11. 

After the data point, x i has been quantized by this method, 
it is replaced by its representative alphanumeric character and 
plotted or printed. The plotting is output to the Stromberg- 
Carlson 4020 recorder for display. 

If automatic quantization is not desired, the program pro- 
vides an option where the user can input quantized levels as a 
table if other methods are desired (reference ERFDP users manual, 
section II of this report) . 

B. Data Problem Parameters 

NCH Number of channels on input tape 

NSPS Number of samples per record or scan 

NSCANS Number of scans to process 

NSKIP Number of scans to skip before processing 

NSTART Starting sample number in the scan 

NSTOP Stopping sample number in the scan 

ITERM Number of passes to make processing 120 samples 
each pass 

N Number of quantized levels plus one 

ICHAN Selected channel used in grey level mapping 
IPRT Option to print grey level 
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IPLT Option to plot grey level 

INCX Increment in x direction for each sample 

(rasters) 

INCY Increment in y direction for each sample 

(rasters) 

NSTX Starting x coordinates or plot frame 

NSTY Starting y coordinates on plot frame 

NBTLG Bit length of input data words 

MODE Signifies FORTRAN or non- FORTRAN 

ITYPE Type of input data; Floating point or fixed 
point 

MSFC MSFC scanner format option 

NCRE Data increment 

LTN Logical tape unit to load input data 

IOPT Option to select automatic quantization 

or input quantized table 

3. CONTOUR PLOTTING 

In the analysis and interpretation of earth resources data, 
a useful tool is means of locating and outlining borders of 
ground scene images and also plot the snow pack temperature pro- 
files recorded by an electronic scanner. A computer program was 
written to generate line plots, which is a line connecting points 
of boundaries, specified altitudes, temperatures, etc., provid- 
ing a graphical display of contour levels or borders of ground 
scene images. 

The data array to be plotted consists of the input data 
array only. These data points are ordered in an x y coordinate 
system. Four adjacent data points and their coordinates are 
examined to determine if a specified value intersects any of 
these four points. If an intersection occurs, these coordinates 
are converted to plotting raster counts, and two points are 
flagged. One point being coordinates of the entry point to the 
four adjacent data points and the other being the coordinates 
of the exit point of the four adjacent data points. If the inter- 
section of a specified level continues to the next four successive 
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adjacent data points, then the coordinates of the exit point 
become the coordinates of the entry point of the next four adja- 
cent data points. A search continues for the coordinates of 
the exit point and this procedure continues until the data set 
is exhausted. Lines are drawn connecting the points, reflecting 
continuity of the specified level. If the entry and exit point 
of successive four adjacent data points exist, then a line con- 
necting these points is continuous. Otherwise, there exists a 
discontinuity and the points will not be connected. Continuous 
lines of multiple levels can be drawn reflecting contours of 
altitude, temperatures, and homogeneous area boundaries. 

A. Computer Program 

The computer program was written as a module to be included 
in the ERFDP and was designed to contour data sets of infinite 
lengths. However, only 2500 data points reside in the computer 
at one time because of the physical storage limits. Since only 
one block of data is processed at one time, the program auto- 
matically reloads blocks of data and abuts each block to provide 
a continuous plot of contours. 

In Figure 2 contours of boundaries were plotted of the first 
120 scans of field Cl from the Purdue data set. The data were 
blocked such that contours of 1776 resolution elements were plotted 
which covered eight scans. Each resolution element covered eight 
rasters on the plot frame, therefore each block abutted together 
occupies 64 rasters. After the abutment process, the boundary 
contours appear continuous. The boundaries that were contoured 
are shown in Figure 3. The boundary contours are elongated slightly 
due to scaling on the plot frame. 


8 



2>0» 


SCAN 1-120 



„C~ v.^~ <r\ • 

^ v yy* 


"S'-r- 

J* / 

<5 5 0 


r^> s’ in 

S' - „ a =xo « vr. ~^2 

jrcri.*-xLy oxsrvia 0 

V „L/ Vn 5 °^ 1 W 

>> . o # t /}m r 

<3 ^ d^f 1 

<=> ?«._> „ tKr: s 
->0 ^XL^A-rC 1 •= 

J? l 

a ^ w A9Z C o 

XV ! °0 

j'.n'HV i 

_/<?o0o V *« 

. / &TJ?y* < r?7 o 

JUSKT V. 

A*;*// ys 2 — c 
s+ V/ 3 /X SJ — ‘-'” 3 
orij' Zr di 

, r-H^ 9 (V ?/ 

%// y • a 5? 

)V JS Cj 

- x sc 


Fig 


C 


RESOLUTION ELEMENT 1-222 




Z >000 


RESOLUTION ELEMENT 1-222 


SCAN 1-120 





RESOLUTION ELEMENT 


Figure 3 


Boundary Determination Channels 4 of 12 


(Purdue Cl) 


B. Data Problem Parameters 


NCH Number of channels on input tape 

NSPS Number of samples in a scan 

IRW Logical tape units to load input data 

NCHAN Channel number 

NSNCRE Number of scans to increment 

NPCRE Number of samples to increment 

NPTSL Lower starting point 

NPTSU Upper stopping point 

NOTLG Bit length of input data word 

ITYPE Type data on input tape; Floating point 
or fixed point 

MSFC MSFC scanner format option 

NSKIP Number of records to skip 

MAXSCN Total number of scans to process 
NSECT Number of sections to process 
MSZX Data block size in x (samples) 

MSZY Data block size in y (scans) 

BLK Plot frame block size in rasters 

FHINC Contour level increment 

ZMIN Minimum contour level 

ZMAX Maximum contour level 

LAB Label interval 

4. ISOMETRIC DISPLAYS 

On occasion it has become necessary to preview the raw data 
and to determine the physical appearance to better understand 
the behavior of data collected from various ground scenes and 
from different types of sensors. For this purpose an isometric 
program was written to produce a two-dimensional projection of 
a three-dimensional graph of the raw data amplitude versus the 
x and y ground scene coordinates. 


11 



A. Computer Program 

Program Isometric will handle any FORTRAN IV formatted 
tape with up to 12 channels of data. Only one channel is plotted 
per computer pass and the channel selection is based upon an in- 
put parameter. 

The program was developed to produce two-dimensional data 
projections of ground cover scenes where data are gathered by 
sensors aboard an aircraft. The data collected across the flight 
line represent resolution elements, and the data collected along 
the flight line represent scans . Any NXM array can be displayed 
where N is the position in an array and M is the value to be dis- 
played. The program does not require large amounts of core 
storage since it performs a continuous operation on the NXM array 
which is refilled before each operation. 

The Stromberg-Carlson 4020 frame reference coordinates are 
altered after each operation for proper scaling based on the input 
parameters provided by the user. These parameters also control 
the density of the points plotted and the angle of rotation de- 
sired for display. Multiple plot frames are generated as neces- 
sary to display all the input data. To provide frame abutment 
for a continuous plot, the last row of each frame and the last 
sample of each row is stored and repeated on the next frame. 

B. Data Problem Parameters 

NCH Total number of channels per resolution 

element on input tape 

NSPS Total number of resolution elements per 

channel across the ground scene 

IRW Logical tape number for input data 

NCHAN Number of the channel that isometrics are 
to be displayed 

NSNCRE Number of scans to increment along the 
flight line; This provides an option to 
use every scan line, every other scan 
line, etc. 
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NPCRE 

NPTSL 

NPTSU 

MAXSON 

YMIN 

YMAX 

NBLSZX 

NBLSZY 

NSECT 

NSMOV 

NDIREC 


Number of resolution elements to increment 
going across the flight line; If every 
resolution element is not desired, resolu- 
tion elements can be skipped. 

Starting resolution element number for 
the isometric 

Stopping resolution element number for 
the isometric 

Total number of scans along the flight 
to display 

The minimum value of the input data used 
for calculating scale factors for plotting 
the y axis 

The maximum value of the input data used 
also for calculating scale factors for 
plotting the y axis 

Stromberg-Carlson 4020 reference frame 
coordinate increment in the x direction 

Stromberg-Carlson 4020 reference frame 
coordinate increment in the y direction in 
raster counts (normal range 6 to 12) 

NBLSZX and NBLSZY are used to determine 
the degree of rotation of the isometric. 

Number of passes through the data neces- 
sary to display the full-scan width; 

If NPTSL and NPTSU only cover a portion 
of the data then NSECT can be adjusted 
to cover all the data. 

Example: 

If NPTSL = 1, NPTSU = 128 and NSPS - 256, 
then NSECT = 2 will cover resolution 
element 1 through 128 and resolution 
element 129 through 256. 

Number of points in a moving mean span 
used in smoothing the input data; Set 
to zero if smoothing is not desired. 

Direction of the rotation of the iso- 
metric; 1 = counterclockwise and 
-1 = clockwise rotation. 
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A portion of flight line Cl (Purdue data set) was used 
in the following examples. Every resolution element and every 
scan was used requiring two passes through the data. NSECT 
was set to 2, but only one section is shown here. NBLSZX and 
NBLSZY were set to 8 and 12 respectively and both counterclock- 
wise and clockwise rotations were used. 

Figure 4 shows scan lines 86-129 and resolution elements 
1-111 of the total 222 elements of channel 8. NDIREC was set 
to -1 to rotate clockwise. Figure 5 shows channel 8 with only 
NDIREC changed to +1 to rotate counterclockwise. Figures 6 
and 7 reflect identical isometrics using channel 2 of the same 
data. 


5. DYNAMIC JOINT PROBABILITY DISTRIBUTION 

Extensive data analysis of data on hand has revealed inade- 
quacies in the current joint probability distribution program. 

The program was restrictive in that a fixed amount of square 
storage array allocation was required. This meant that, the 
greater the spread of the clusters of the joint probabilities 
or individual data pairs, the larger the square storage array 
required. This also required a screening of all the data to 
determine the ranges of data pairs. After determining the data 
ranges, the minimum square storage array required to display 
the joint probabilities became NXM, where N is the range between 
the minimum and maximum value of the X data in the joint pairs, 
and M is the range between the minimum and maximum value of the 
Y data in the joint pairs. As this clearly points out, this does 
not adapt itself favorably to data that is ill behaved and pos- 
sessing widely spread data clusters or widely spread individual 
data pairs. Due to the physical limits of the computer storage, 
in numerous cases, a considerably amount of data was lost. This 
depended largely on the characteristics of the data and since 
the primary function of the program was to display the dependent 
characteristics of two data channels, this required rerunning 
the program and making adjustments to the limits. Therefore, 
the program became inadequate. 
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Figure 5 Flight Data Display 
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Figure 6 Flight Data Display 
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Figure 7 Flight Data Display 



A new technique has been employed, that adapts itself to 
data of any nature, particularly cluster data where the data 
clusters are tremendously scattered. The technique requires 
only as much storage as actually required by the number of dif- 
ferent data pairs and the size of the data storage required is 
not dependent on the ranges of the data. 

6. JOINT PROBABILITY DENSITY FUNCTION 


A joint probability density function program is a necessary 
tool in the analysis of earth resources flight data. This pro- 
gram gives a preliminary insight to the commonality of occur- 
rences observed from two time-series data traces. 

For applications to feature extraction, the outputs of this 
program give insight to the distinctness of different features 
and the location of decision boundaries necessary to separate 
different features. 

The joint probability distribution is calculated digitally 
by selecting a range [a, b] on both sets of data X, Y and divid- 
ing the range into K intervals. This gives 


C 


b-a 

K 


= scale factor 


To find the X and Y coordinate for each variable X and Y jointly, 
an indexing pointer is calculated for each value of X and Y by 


I 


x 


(i V a > 

c 


and I y 


(r n - a > 

c 


An integer one is added to the coordinate P(I x »I y ) for every value 
of X and Y in the interval [a,b]. The coordinate P(I x >Iy) con- 
tains the total number of occurrences. 
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A. Computer Program 


Data points are paired (user option), given a class, and 
are stored in two single arrays, one for the X coordinates, 

NP X , and one for the Y coordinates, NP^. The number of data 
pairs that are common is accumulated in an array NKNT (i) iden- 
tifying the class and the number of occurrences, i, of that class. 

The first data pair read from the input tape, will be as- 
signed class number 1 automatically. All subsequent data pairs 
read from the input tape, will be compared with all the indivi- 
dual data pairs that have been assigned a particular class. If 
no match is found during the comparison, a new class is created, 
using this data pair. There can be up to 4000 different classes 
of paired data points stored. 

Once all the data pairs are classified, the data pairs are 
arranged and sorted such that the Y coordinates or vertical data 
are sorted in descending order and the X coordinates or horiz- 
ontal data are sorted in ascending order. Data are rearranged 
in such a way that the largest spread between the minimum and 
maximum is displayed on the vertical axis, which is printed down 
the print paper. An example of the three data arrays containing 
the joint probability before and after the vertical sort is shown 
in Figure 8a and 8b. 


B. Core Storage Image 

For any data configuration, the core storage image will 
appear as in Figure 9a. Projecting this onto an NXM storage 
array used in the conventional method, shows the wasted storage 
which is used to store blanks (or no occurrences). See Figure 9b. 
The conventional data storage method required for this example 
is 6x7, or 42 core storage locations. The present scheme re- 
quires only 21 core storage locations to display the above example. 
The advantage, also, is that only 21 core storage locations are 
necessary regardless of the scatter of the joint probabilities. 
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In Figure 10a and 10b is an example of the same number of paired 
classes, but the data pairs are more widely scattered. In this 
example, the conventional data storage required is 14 x 12 or 
168 core locations. The present scheme still required only 21 
core locations. 


C. Output Display 

To display the arrays as a joint probability distribution, 

the first location in the NP^ array is examined and a decremental 

counter is set equal to this value. All the data points in the 

NP array that are paired with this NP value, are collected 
x y 

along with their number of occurrences NKNT(i), and stored in a 
working array. This working array is sorted in ascending order 
for printing from left to right across the print paper. A print 
line is loaded with blank characters and the blank character is 
replaced with the number of occurrences in that location, if any 
exist. The number of occurrences is designated by an alphanumeric 
character of some hierarchy ordered by the user. The print line 
is then output to the printer. The next location in the NP y 
array is examined and the counter is decreased by one. If the 
NPy value is less than the counter, a print line filled with 
blank characters is printed and tagged with the value of the 
counter. The counter is then decremented and compared with the 
NPy value again. If there is a comparison, then the above proce- 
dure is repeated. This continues until the NP^ array is exhausted. 

D. Data Problem Parameters 

NCH Total number of channels on input tape 

NSPS Samples per logical record or scan line 

NS CANS Number of logical records or scan lines 

to process 

NSKIP Initial physical records to skip before 

processing 

NSTART Starting sample number 
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NSTOP 

Stopping sample number 

LTN 

Logical tape unit to load input data 

NOJP 

Number of paired joint probabilities 
to process 

IMX 

X channel select 

IMY 

Y channel select 

SCALE 

Used to scale data 

BIAS 

Used to shift data 


7. BOUNDARY MAPPING 

A concerted effort has been made to improve boundary mapping 
techniques and to extend these ideas and techniques in order to 
classify different homogeneous areas of ground scene data. This 
section briefly describes one of several boundary enhancement 
techniques, which were investigated, and contains examples of 
test cases and problem parameter inputs. 

The technique employed in this computer program, incorporates 
a moving rectangle made up of four adjacent resolution elements 
moved successively through the data. The configuration, produced 
by the magnitude of any four adjacent resolution elements con- 
nected by an imaginary path, will be defined as the area. Iso- 
metric displays of the ground scene (Figure 11) show different 
area configurations produced by any four adjacent resolution 
elements. The possible area configuration models are shown in 
Figure 12. From these configurations, the equations for calculat- 
ing the area were derived, based on two cases presented by the 
data. Case I is the area produced by four adjacent resolution 
elements along the Y and Z plane, and case II is the area pro- 
duced by four adjacent resolution elements along the X and Z 
plane. The equations that were derived (A-^, B-^, for case I 
and A£, ^2 ^ or case H) are shown in Figure 13a. 

Several equations that appear here are redundant and are 
eliminated by combining similar equations. The composite test 
statements that determine which equation to use, are shown in 
Figure 13b, and a flow diagram of the decision logic incorporated 
in the computer program, is shown as Figure 14. 
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Figure 12 Area Configurations Produced by Four Adjacent Data Samples 
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Figure 13 Area Equations for Four Adjacent Data Samples 
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Figure 14 Flow Diagram of Decision Logic 



A probability density distribution of the area is calculated 
for each channel of data and every scan. An estimate of this 
probability density function is obtained digitally by dividing 
the range for X into an appropriate number of class intervals, 
say d^, and by tabulating the number of occurrences in each of 
these intervals. The number of occurrences in each interval 
satisfies the equation 

k-1 

Ni 

i=0 


N, 


N, 


of 

class 

intervals 

and 

N is 

the 

sequence 

{ N i> 

is 

found by 


= 

Total 

number 

of 

X: 

x C 

d l 

= 

Total 

number 

of 

X: 

d l< 


=5 

Total 

number 

of 

X: 

d 2 < 



N k = Total number of X: d^^ < X <£_ d^ 


A plot of the probability distribution of channel 3 and 
channel 12, both horizontal and vertical, are shown in Figures 15 
and 16. These plots reflect the distribution of one scan of data 
in particular, scan number 5. The criterion to determine if a 
resolution element reflects a change in a homogeneous population 
of ground scene is the mean of the probability distribution of 
the area plus its associated 6. All area calculations of four 
adjacent resolution elements that fall outside of the mean plus 
6 will be flagged as a boundary. 
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Figure 15 Boundary Determination, Probability Distribution Channel 3 



Figure 16 Boundary Determination, 


Ch 12 



Distribution Channel 12 


The mean area of the probability distribution, P(A) , is 
calculated by 


PW 



i=l,2,.. .MAX 


on 6 by 



** 1/2 



i-1,2. . .MAX 


In the test cases, all 12 channels were used, except for 
two cases where channel 3 and channel 12 were run separately, 
and boundary flags set for each channel. The final decision 
for determining the boundary, is based on an input parameter 
selected by the user. The number of channels indicating a 
boundary has to be greater than this value. In Figure 17 only 
channel 12 was used. This produced too many boundaries and 
the separation of different homogeneous areas was not clearly 
defined. In Figure 18 only channel 3 was used; this produced 
too few boundaries. Figure 19 shows all 12 channels with at 
least seven reflecting boundaries (Input Option); this produced 
too few boundaries. Figure 20 shows at least three channels of 
the 12 reflecting boundaries which produced too many boundaries. 
Figure 21 produces the best boundary map which uses at least 
four channels of the 12 total. 
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Figure 19 Boundary Determination Channels 7 of 12 (Purdue Cl) 
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A. Data Problem Parameters 

Number of scans or logical records to process 
an input tape 

Starting resolution element in the scan 

The number of resolution elements or columns 
to process from each scan (<^ 255) 

Number of channels present on input tape 

Number of variables or intervals desired in 
calculating the probability distribution; 
Maximum of 2000. 

Number of alphanumeric characters to display 
boundaries (note: At the present only two 

characters are used, blanks and dots, no 
boundary or yes boundary respectively) . 

Total number of channels used in initial 
boundary calculation 

Number of channels (plus one) , of ISUM, 
necessary for a final boundary decision 

Number of bits per word on the input tape 

Mode of input tape (2=non- FORTRAN, 1=F0RTRAN) 

Number of records to skip before processing 
input data 

Data increment 

Number of raster counts per character desired 
in displaying boundary on the Stromberg- 
Carlson 4020 plotter recorder 

Raster counts incrementation in the X direction 

Same as above in the Y direction 

Starting position or frame along the X axis 

Starting position on frame along the Y axis 

Channel number desired in calculation of 
boundary 



8. CHANNEL ALIGNMENT 


In recent months there has been an increasing demand to 
process and analyze data collected by array cameras using dif- 
ferent filters in the frequency spectrum. In order to analyze 
data from this camera source, the camera film is digitized and 
the digitizing process is carried out separately for each fil- 
tered film. Gross errors are introduced into the data, mainly 
missalignment of the film when being digitized by a densitometer. 
This missalignment must be compensated for before the digitized 
data can be merged and processed. A process was developed to 
align the digitized data sets by a matching matrix technique. 

Boundaries separating homogeneous areas are calculated from 
each data set (see section 7) and the boundaries of each data 
set are correlated and matched. One particular technique pro- 
posed was the use of NXM data matrix for each data set where one 
set is assigned the reference channel and the other sets are 
correlated with the reference channel by moving the matching 
matrix successively through the selected data area. 

The correlation of the reference channel with other channels 
is done by matching or correlating a resolution element from the 
reference channel with all resolution elements in the matching 
matrix of other channels. If there exists a match or correlation 
of boundaries, then a cell in the matching matrix display at those 
coordinates is incremented. The NXM matching matrix is moved by 
one resolution element, and the next resolution element from the 
reference channel is matched or correlated with all the resolution 
elements in the matching matrix of other channels. If there exists 
a match or correlation of boundaries, then a cell in the matching 
matrix display at those coordinates is again incremented. This 
continues as the matching matrix is moved successively through 
the selected data interval. 
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An expanding and collapsing technique is employed at the 
beginning and end of the data set to eliminate edge effects and 
loss of data. As the NXM matching matrix is moved successively 
through the data set, the correlation of the two channels, if a 
correlation exists, is accumulated at the beginning of the data 
set at the matching matrix display coordinates MATRIX (i,j) 
going across, where 

i = | - b+1, .N - b+2, . ..N; k = 1,2,3. . , 

and going down where 

j * ^ ■ b+1, - b+2 , . . .M; k — 1,2.3...?^ . 

This continues until the entire matching matrix is within the 
objective data, thus eliminating the edge effect. Once the match- 
ing matrix is within the objective data then i = 1,2,3...N; and 
j = 1,2,3...M. Upon termination of the matching matrix at the 
end of a data set, a collapsing technique is employed. This tech- 
nique is the same as the expanding technique except the accumula- 
tion at the end of the data set at the matching matrix display 
coordinates MATRIX (j,i) is i = 1,2.3... N-k, where k = 1,2,3... 
and j = 1,2,3... M-k, where k = 1,2.3... Ir • 

The matching matrix is then printed for each combination of 
correlations. No boundary and boundary correlation of two data 
channels (0.1), and boundary and boundary correlation of two data 
channels (1,1). Alignment shifts can be performed by locating 
the peaks in the NXM matrix where boundaries of two channels are 
matched. The matching matrix is shown in the output example. 

This peak occurred at coordinates (-5, +4), therefore, the channel 
being correlated with the reference should be shifted -5 data 
records and +4 resolution elements to align the two channels. 
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This technique corrects for any horizontal and vertical 
alignment satisfactorily, however it does not correct for any 
skewness introduced, while the film is being digitized, or any 
distortion in the camera array angles. A technique that fully 
corrects for this type of roissalignment is presently being 
developed, but due to inadequate computer turnaround time, re- 
sults are not available. 

A. Computer Program 

A program was written to match and correlate up to three 
data channels and as many different intervals of a data set as 
desired. This will enable the user to calculate the alignment 
in areas of the data set where boundaries are clear and can be 
detected. Input to the program is three FORTRAN formatted fixed 
point binary tapes containing boundary flags (fixed point integer 1 
for boundary and 0 for areas that are homogeneous) . 

B. Data Problem Parameters 

NTIMES Number of sections to correlate 

LT9 Logical unit to load input tape number one 

LT10 Logical unit to load input tape number two 

LT11 Logical unit to load input tape number three 

NSTART Starting resolution element 

NSTOP Final resolution element 

NSCANS Number of scans in the matching matrix 

NSPS Number of resolution elements across the scan 

in the matching matrix 

NREFX Starting resolution element in the reference 

channel; Tape number one. 

NREFY Starting scan number in the reference channel; 

Tape number one. 

NCHlX Starting resolution element in the channel to 
be correlated; Tape number two. 

NCHlY Starting scan in the channel to be correlated; 

Tape number two. 

NCH2X Starting resolution element in the channel to 
be correlated; Tape number three. 

NCH2Y Starting scan in the channel to be correlated; 

Tape number three. 
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C. Program Flow Chart 
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Channel Alignment (Concluded) 
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9. CHANNEL REGISTRATION 

Registration of digital images of the digitized film pro- 
cessed from a multiple camera array system, has become necessary 
due to the rotational, translation and scaling errors introduced 
during processing. A course alignment is made by employing the 
matching matrix technique but this only corrects for a horizontal 
or vertical shift in the digital images for alignment. There 
exist errors in the digital images, such as rotation and transla- 
tion errors, which have to be corrected by a more effective means. 

A technique which will accomplish this is presently being developed 
and has been demonstrated. Computer results are not available, 
due to contract expiration which prevented finalizing the analysis, 
therefore, no illustrations are shown. 

The technique requires generating boundary maps of the digital 
images and a manual technique is used to correlate similar areas 
of the digital images with each other. A good sampling of the 
entire scene is made which provides enough input to a set of equa- 
tions (equations are presently being derived and improved upon) 
to calculate coefficients for each digital image combination. 

When used with the x and y coordinates of a reference channel, 
position pointers are calculated giving the respective x and y 
coordinates of the channels to be aligned. Data points are taken 
from these coordinates and placed at the coordinates of the ref- 
erence channel. All channels are then output to magnetic tape 
becoming merged for further processing. 

A program was written to overlay any combination of digital 
images containing boundary information (0 and 1) that have been 
registered and merged. This overlay program produces a map showing 
the location of boundaries and homogeneous areas that are similar. 
This program is used in conjunction with the registration process 
and evaluating the accuracy of the registration process. 

After the digital images are aligned satisfactorily the reg- 
istration process can be carried out on the raw data and processed 
as desired, or the merged boundary tapes can be processed. 
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A. Computer Program 


A computer program was written to register multiple digital 
images using boundary information or raw data. The program inputs 
a reference channel and the x and y coordinates of data points 
from the reference channel are input to a subroutine FETCH which 
returns through the call statement the x and y coordinates of the 
channel to be aligned. Subroutine FETCH contains DATA statements 
with coefficients used in the calculations. These DATA statements 
can be changed to reflect coefficients of different digital images. 

B. Data Problem Parameters 

NSPS Number of samples per scan 

NCHAN Number of channels on input tape 

NSCANS Number of scans to process 
NSKIP1 Initial position of tape number one 
NSKIP2 Same as above for tape number two 

NSKIP3 Same as above for tape number three 
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C. Program Flow Chart 



Channel Registration 
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Channel Registration 
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10. SPECTRAL DISCRIMINATION 


Capabilities have been developed to discriminate ground 
scene features using computer techniques with digitized images. 
These capabilities are available at MSFC for analysis and inter- 
pretation of earth resources flight data. Computer programs and 
algorithms were developed during this contract period, March 22, 
1971 to September 22, 1972, under the auspices of the National 
Aeronautics and Space Administration, with Mr. Robert Jayroe 
COR. The mathematical rationale used in the development of these 
computer programs are to be published by Mr. Robert Jayroe in 
a NASA report entitled "Unsupervised Spatial Clustering with 
Spectral Discrimination." 
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ERF DP USERS MANUAL 
SECTION II INTRODUCTION 

Processing earth resources data is the future will require 
a somewhat automatic approach due to the volume of data and the 
extent of analysis required to interpret the data. This requires 
a complex data processor designed to maintain an automatic data 
processing (ADP) environment for the user. The processor is 
designed using independent modules with an executive program 
driving each module. A block diagram depicting this modular 
structure is shown as Figure 22. 

This processor uses an overlay technique, which compiles 
and communicates with all modules, but only executes modules 
selected by the user . The concept of using overlays provides 
capabilities for processing input data, through all the program 
modules sequentially or only selected modules as desired. Each 
module can be improved upon, expanded and modified independently 
of other modules , and each module can be improved or modified 
in a 32K core storage environment irrelevant to the other modules 
(see Figure 23) . 

This conept has the disadvantage of having to utilize an 
intermediate storage tape, which necessitates executing read/write 
commands, which are certainly time-consuming. However, the ad- 
vantages outweigh the disadvantages significantly by providing 
versatility in the overall function. To further improve on the 
e ffi- c i enc y the program, the read/write commands use alternate 
input/output channels when possible to eliminate, in most cases, 
the computer being in a "wait" condition due to either channel "A" 
or channel "B" being busy. 

The overlay concept entails having a main program which calls 
each individual module as requested by the user. This keeps the 
overall processing contained in one large program and each program 
module is brought into the computer, compiled, overlay structures 
created, and then output onto an overlay tape. All modules are 
treated in this fashion regardless of which and how many are to 
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Figure 22 Block Diagram Depicting Modular Structure 



















be used. Each module to be executed is brought into the computer 
in a "transit" area (see Figure 24), and upon completion, another 
module, if requested, is brought into the same "transit" area to 
be executed. A common area in the computer is used to store in- 
formation that may be used by more than one module. Each module 
has access to this common area at any time it occupies the 
"transit" area. 

1. EARTH RESOURCES FLIGHT DATA PROCESSOR 

The ERF DP is comprised of seven separate and independent 
modules which reside in computer memory only while being executed. 
Module seven consists of three lower level modules that are called 
into memory to perform their function automatically when module 
seven is requested by the user. The segmented structure appears 
as diagrammed. 



a 



b 


b 


Each module can be called into core memory and executed 
independently and in any order with the exception of module seven. 
Module six has to be executed prior to module seven, since the 
boundary mapping output is input to module seven. If module six 
has been run previously and a boundary mapping output tape is 
available, then module seven can be executed alone. 


56 




TAPE 


PRINT 


CARDS 


GRAPH 


Figure 24 Block Diagram ERFDP 

Earth Resources flight Data Processing 


57 









A dimensional variable called MODULE is dimensioned 8, 
which sets aside eight memory locations to store a users request 
for module execution. The number of locations filled depends 
upon the number of modules requested to be executed by the user. 
The remaining locations contain zeroes and when the processor 
encounters a zero while examining the list, it terminates under 
system control . 


A. Data Problem Parameters Example Setup 


$DATA 

$INIT 

MODULE = 1,2, 4, 6, 7 
$END 


User request modules 1,2,4, 6, 7 
to be executed and in that order. 


B. Input Tapes 

Unit - Users option under individual modules (A6 recommended) 
Type - Reference individual modules 

C. Output Tapes 
Reference individual modules 
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2. MODULE ONE 


Module one when called into memory calculates a probability 
density function for any number of channels not to exceed 12. 

The module provides options to calculate probability density 
functions on selected channels also. 


A. 


$ INPUT 1 
NCH = 3 
NSPS = 794 
N SCANS =500 
NSTART = 1 
NSTOP =794 
NBTLG = 12 
MODE = 2 
I TYPE = 0 


MSFC = 0 

LTN = 10 
NSKIP = 2 
NCRE = 1 
XMAX = 1200.0 
XMIN = 300.0 
NOCHS = 3 


NWHICH = 1,2,3 

$END 

7/8 


Data Problem Parameters Example Setup 


Number of input channels 

Number of samples per record 

Number of scans or records to process 

Starting sample number 

Stopping sample number 

Bit length of input data word 

Signifying non- FORTRAN input tape 

Fixed point input; ITYPE = 1 ; Floating 
point 36-bit word input 

Not MSFC scanner format; MSFC = 1; MSFC 
scanner format 

Logical tape unit of input tape 
Skip first two data records 
Data increment to process 
Maximum value in data set 
Minimum value in data set 

Number of channels to calculate probability 
density function 

Calculate probability density function on 
channels 1,2, and 3 


B. Input Tapes 

Unit - Users option under input parameters 

Type - Any odd parity binary, 3-bit modulus, fixed point 
with word lengths 36 bits, or floating point. 


C. Output Tapes 
None (only print output) 
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D. Program Flow Chart 
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Multiple Channel Probability 
Distribution (concluded) 



OUTPUT EXAMPLE 

PROBABI LITY DISTRIBUTION 


AMPLITUDE 

CHI 

CH 2 

CH 3 

CH 4 

CH 5 

CH_6 

CH 7 

CH 8 

CH 9 

CH 10 

CH 11 

CH 12 

-1.5 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1.4 

0 

0 

2 

0 

2 

0 

0 

0 

0 

0 

1 

0 

-1.3 

0 

0 

3 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1.2 

1 

0 

9 

0 

5 

0 

0 

1 

0 

0 

1 

0 

- 1 .] 

0 

0 

21 

0 

7 

0 

0 

0 

1 

1 

8 

0 

-1.0 

5 

0 

29 

0 

15 

1 

0 

1 

0 

0 

22 

0 

- .9 

4 

0 

35 

0 

17 

1 ■ 

0 

3 

2 

4 

31 

0 

- .8 

8 

1 

39 

0 

34 

2 

3 

9 

7 

8 

47 

0 

- .7 

12 

0 

48 

0 

39 

8 

8 

18 

20 

16 

60 

0 

- . 6 

24 

7 

55 

1 

47 

12 

12 

22 

25 

30 

55 

0 

- .5 

21 

11 

66 

0 

51 

22 

22 

38 

40 

40 

40 

0 

- .4 

38 

5 

58 

2 

55 

31 

41 

47 

48 

52 

27 

0 

- .3 

• 45 

12 

49 

12 

49 

49 

70 

50 

55 

60 

11 

0 

- .2 

57 

22 

42 

18 

44 

66 

82 

54 

62 

72 

6 

5 

- .1 

59 

37 

33 

20 

39 

72 

96 

47 

51 

55 

5 

1 

0 

62 

44 

25 

35 

33 

62 

79 

41 

45 

52 

0 

8 

.1 

58 

48 

19 

30 

25 

40 

62 

35 

39 

31 

1 

16 

.2 

51 

49 

12 

41 

18 

22 

38 

25 

25 

20 

0 

27 

.3 

47 

56 

6 

54 

7 

12 

12 .. 

17 

16 

15 

0 

32 

.4 

35 

52 

1 

61 

8 

8 

10 

12 

8 

7 

0 

39 

.5 

25 

40 

0 

55 

2 

2 

7 

8 

5 

3 

0 

47 

.6 

22 

32 

1 

45 

1 

0 

2 

0 

1 

5 

0 

52 

.7 

8 

24 

0 

40 

0 

1 

1 

2 

0 

1 

0 

65 

.8 

t „ 

16 

0 

35 

1 

0 

0 

0 

0 

0 

0 

72 

.9 


8 

0 

28 

0 

0 

0 

0 

0 

0 

0 

58 

1.0 


7 

0 

20 

0 

0 

0 

0 

0 • 

0 

0 

46 

1.1 

1 

1 

0 

12 

0 

0 

0 

0 

0 

0 

0 

39 ' 

1.2 

0 

0 

0 

7 

0 

0 

0 

0 

0 

0 

0 

22 

1.3 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

7 

1.4 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

1.5 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


NUMBER OF OCCURRENCES 



3. MODULE TWO 


Module two outputs on print paper and/or Stromberg-Carlson 
4020 recorder, the quantized levels of a data set with alpha- 
numeric characters displaying specific quantized levels. Alpha- 
numeric characters can be selected and input to the module on a 
card included in the input problem parameters to provide different 
shading for certain levels. This module displays only one channel 
per computer pass, and that one channel is optional and selected 
by the user. 


A. Data Problem Parameters Example Setup 


$INPUT2 
NCH = 3 
NSPS = 794 
NSCANS = 500 
NSKIP - 2 
NSTART = 1 
NSTOP = 120 
ITERM = 6 

N = 11 
ICHAN = 3 
IPRT = 1 
IPLT = 1 
INCX = 0 

1NCY = 8 

NSTX = 0 
NSTY - 0 
NBTLG = 12 
MODE = 2 
ITYPE = 0 

MSFC = 0 


Number of channels on input tape 

Number of samples per record or scan 

Number of scans to process 

Number of scans to skip before processing 

Starting sample number 

Stopping sample number 

Number of passes processing 120 samples 
each pass 

Number of levels plus 1 
Channel used in mapping 
Option to print 
Option to plot 

Increment in x direction for each sample 
(rasters) 

Increment in y direction for each sample 
(rasters) 

Starting x coordinates on plot frame 
Starting y coordinates on plot frame 
Bit length of input data words 
Signifies non- FORTRAN 

Fixed point input data; ITYPE = 1 ; Input 
data floating point 36-bit words 

Not MSFC scanner format; MSFC = 1; MSFC 
scanner format 
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NCRE = 1 
LTN = 10 
IOPT - 0 
$END 
$GRYLYL 
GLVL = 1. 

1 . 


1 . 

$END 

7/8 


Unit 

Type 


None 


Data increment 

Logical tape unit of input tape 
Signifies automatic quantization 



o 


o 


2., 3., 4., 5., 6., 7., 8., 9., 10., 11. 


B. Input Tapes 

- Users option under input parameters 

- Any odd parity binary, 3- bit modulus, fixed point, 
with word length <£_ 36 bits or floating point. 


C. Output Tapes 
(only print output) 
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D. Program Flow Chart 
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Grey Level Mapping 
(cone luded) 




OUTPUT EXAMPLE 
PURDUE FLIGHT LINE Cl 
AUTOMATIC QUANTIZED GREY LEVEL 
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OUTPUT EXAMPLE 
PURDUE FLIGHT LINE Cl 
AUTOMATIC QUANTIZED GREY LEVEL 


RESOLUTION ELEMENT 112-222 SCANS 1-120 
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4. MODULE THREE 


Module three calculates contour lines throughout an input 
data set depicting elevation, altitude, temperature or boundary 
contours. The contour lines are output and displayed by the 
Stromberg-Carlson 4020 plotter. Frame butting in the module pr 
vides for continuous contour plotting for infinite data sets. 
The contour lines can be scaled by utilizing the block size 
option (BLK) in the input problem parameters. 


A. Data Problem Parameters 


$INPUT3 


NCH = 1 

Number of channels 

NSPS - 90 

Number of samples in a scan (record) 

IRW = 11 

Logical tape unit for input data 

NCHAN = 1 

Channel number 

NSNCRE = 1 

Number of scans to increment 

NPCRE = 1 

Number of samples to increment 

NPTSL = 1 

Lower starting point 

NPTSU = 90 

Upper stopping point 

NBTLG =36 

Bit length of input data word 

MODE = 1 

FORTRAN formated (1) or non-FORTRAN (2) 

ITYPE = 0 

Fixed point input ITYPE = 1, 


floating point input 

MSFC - 0 

MSFC = 1; MSFC scanner format 

NSKIP = 0 

Number of records to skip 

MAXSCN = 100 

Total number of scans to process 

MSZX = 90 

Data block size in x (samples) 

MSZY = 8 

Data block size in y (scans) 

BLK = 8.0 

Plot frame block size in rasters 

FHINC =10.0 

Labeling increment 

ZMIN =8.0 

Minimum label 

ZMAX = 138.0 

Maximum label 

LAB = 1 

Label every interval 

$END 


7/8 
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B. Input Tapes 
Units - Optional under input data 

Type - Any odd parity binary 3-bit modulus, fixed point, 
with word lengths ^ 36 bits, or floating point. 

C. Output Tapes 

SC 4020 S troraberg-Carlson formatted tape 
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D* Program Flow Chart 
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Contour Level Plotting 
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OUTPUT EXAMPLE 


SCAN 1-120 


RESOLUTION ELEMENT 1-222 
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5. MODULE FOUR 


Module four displays a two-dimensional plot of three- 
dimensional data where only one single array of data is input 
and the data samples are ordered in an XY coordinate system, 
where Y is the scan line count and X is the resolution element 
in the scan. Data sets of infinite length can be displayed on 
the SC4020 recorder. The module provides for multiple passes 
through the data set in order to view the entire data set in 
sections . 


A. Data Problem Parameters Example Input 


$INPUT4 
NCH = 12 
NSPS - 222 
NSKIP - 0 
NBTLG - 12 
MODE = 1 

IRW = 10 
NCHAN = 6 
NSNCRE - 1 
NPCRE = 1 
ITYPE - 1 

MSFC = 0 

NPTSL = 1 
NPTSU =111 
MAXSCN = 120 
XMIN =0.0 
XMAX = 222.0 
YMIN = 0.0 
YMAX = 255.0 
NBLSZX 


Number of channels on input data tape 

Number of resolution elements per scan 

Initial records to skip 

Bit length of the input data word 

FORTRAN formatted (1) or non-FORTRAN 
formatted 

Logical tape unit to load input 

Channel selected to be plotted 

Scan incrementation 

Resolution element incrementation 

Fixed point input (0) or Floating point 
input (1) 

MSFC scanner format (1) or not MSFC 
scanner format 

Lower resolution element to start 

Upper resolution element to stop 

Scan lines to be processed 

Minimum value for scaling the X axis 

Maximum value for scaling the X axis 

Minimum value for scaling the Y axis 

Maximum value for scaling the Y axis 

The distance separating the plotted points 
in the X direction 
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NBLSZY 

NSECT = 2 

NSMOV - 0 
NDIREC = 1 

$END 

7/8 


The distance separating the plotted points 
in the Y direction 

Two passes through the data to plot the 
data set in two sections 

No smoothing on input data 

Orientation of the isometric left view (-1) 
or right view (1) 


B. Input Tapes 

Unit - Users option under input parameters 

Type - Any odd parity binary, 3-bit modulus, fixed point 
with word lengths 36 bits, or floating point 
binary. 

C. Output Tapes 

Unit - A8 SC4020 
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D. Program Flow Chart 



Isometric Program Display 
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Isometric Program Display 
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Isometric Program Display 
(concluded) 
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_ OUTPUT EXAMPLE 

FLIGHT LINE Cl 
SCANS 1-43 * 

RESOLUTION ELEMENTS 1-111 
ROTATED CLOCKWISE 
CHANNEL 8 0.58-0.62 (MICRONS) 



RESOLUTION ELEMENT 
DIRECTION 



OUTPUT EXAMPLE 

FLIGHT LINE Cl 
SfANS 1-43 

RESOLUTION ELEMENT 1-111 
ROTATED COUNTER-CLOCKWISE 
CHANNEL 8 0.58-0.62 (MICRONS) 



RESOLUTION ELEMENT 
DIRECTION 



6. MODULE FIVE 


Module five calculates a joint probability density function 
from two time-series data traces selected by the user in the 
data problem parameters. The program provides for multiple passes 
through the data set for multiple joint probability density func- 
tions. The raised data points are classified and stored in one 
single array along with its associated occurrences of like data 
pairs. The joint pairs are sorted vertically in descending order 
and horizontally for each print line in ascending order. 


A. Data Problem Parameters Example Setup 


$INPUT5 
NCH = 12 
NSPS = 222 
NSCANS = 920 

NSKIP = 0 

NS TART = 2 
NSTOP = 222 
LTN = 10 
NOJP = 2 

IMX = 4,3 
IMY * 6,7 
SCALE =1.0 
BIAS - 0 
$END 

Qt ® ),/* ABC... 
7/8 


Total number of channels on input tape 

Samples per logical record or scan line 

Number of logical records or scan lines 
to process 

Initial physical records to skip before 
processing 

Starting sample number 

Stopping sample number 

Logical FORTRAN IV tape unit 

Number of paired joint probabilities to 
process 

Use channel 4 and 3 in X 
Use channel 6 and 7 in Y 
Used to scale 1 
Used to shift data 

(Alphanumeric characters) 
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B. Input Tapes 
Unit - A6 (optional) 

Type - Any odd parity binary 3- bit modulus, fixed point 
with word lengths <£_ 36 bits, or floating point 
binary. 


C. Output Tapes 
(none print output only) 
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D. Program Flow Chart 
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loop for 
no of joint 
distributions 
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input 
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output 
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Dynamic Joint Probability 
Distribution (concluded) 
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7. MODULE SIX 


Module six calculates boundaries outlining edges of homo- 
geneous areas of ground scene data. These boundaries are dis- 
played on print paper where a selected alphanumeric character 
signifies a boundary and blanks signify areas of homogeneity. 

This output is also displayed on the SC4020 plotter where bounda- 
ries are flagged by an alphanumeric "period.” The boundary out- 
put is used in further processing to extract information for 
feature discrimination. 


A. Data Problem Parameters Example Setup 


$ INPUT 6 
NSCANS = 120 
NSTART = 1 
NSPS = 222 
NCH = 12 
NSYM =49 
ISUM = 5 

NBTLG - 12 
MODE = 1 

I TYPE = 0 

MSFC = 0 

NSKIP = 0 
INCX - 0 

INCY = 4 

NSTX = 0 
NSTY = 0 
NCRE = 2 


Process 120 scans of the data set 

Starting sample number 

Number of samples to process; ( 255) 

Number of channels on input tape 

Dimension of the Alphanumeric array 

Number of channels used in calculating 
boundaries 

Bit length of input data word 

FORTRAN formatted (1), non-FORTRAN 
formatted - 2 

Input data fixed point. Floating point 
binary = 1 

Not MSFC scanner format MSFC scanner 
format = 1 

Skip no records before processing 

Incrementation in the X direction on the 
plot frame between resolution elements 
for each scan 

Incrementation in the Y direction on the 
plot frame between resolution elements 
for each scan 

Starting X coordinates on the plot frame 

Starting Y coordinates on the plot frame 

Data incrementation; Use every other data 
sample for boundary calculations 
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$END 

$NCHUSE 

NWHICH = 1,3,4,8,12 Channel selection for calculating 

boundaries 

$END 

©. © *(- =, $+ABC Z12 ...90- 

7/8 


B. Input Tapes 

Unit - Any odd parity binary 3-bit modulus, fixed point 
with word lengths <£_ 36 bits, or floating point 
binary. 


C. Output Tapes 


Unit - B6 

Type - Fixed point binary l's and 0’s FORTRAN formatted 

255 words or less plus 1 FORTRAN index per logical 
record. 
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D. Program Flow Chart 



INPUT 
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Boundary Mapping 
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tou<z 


OUTPUT EXAMPLE 


SCAN 1-120 


RESOLUTION ELEMENT 1-222 



RESOLUTION ELEMENT 




Boundary Determination Channels 4 of 12 (Purdue Cl) 



8. MODULE SEVEN 

Module seven consists of three submodules called into core 
memory sequentially to classify land use data. These modules 
perform spatial clustering, sequential merging, and spectral dis- 
crimination of ground scene images. The module provides an op- 
tion to select multiple passes to further classify small insigni- 
ficant homogeneous areas that may be overlooked or misc lassif ied 
during the first pass. 

The input to module seven is the boundary mapped information 
(see Boundary Mapping) and the raw spectral digitized images. 

The spatial clustering submodule searches areas on the boundary 
map tape that contain homogeneous areas consisting of at least 
100 square resolution elements as a threshold. When this search 
is satisfied, this homogeneous area is identified with a number 
which is incremented for each separate homogeneous area. A new 
tape is created containing these identified clusters. The sequen- 
tial merging submodule searches the tape containing the identified 
clusters and selects the raw spectral data associated with each 
cluster. Criterion parameters are calculated and the merging 
process takes place merging all similar clusters. 

The spectral discrimination submodule, classifies the merged 
clusters starting with an initial class of one, with unit incre- 
mentation for each new class detected. The raw spectral data 
criterion parameters, for each resolution element are compared 
with all classification criterion parameters. If a comparison 
exists, then the resolution element becomes classified. If there 
is no comparison, the resolution element is assigned a blank or 
a boundary flag, (if the resolution element reflects a boundary). 

If the ground scene is not classified satisfactorily, another 
pass can be performed reducing the homogeneous square area thres- 
hold to 36 resolution elements. This will enable smaller homo- 
geneous areas, that were not detected by the 100-square resolution 
element threshold to be detected and an attempt will be made to 
classify these areas also. 
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In the output example of module seven, f'igure 25 shows the 
initial clustering of homogeneous areas that contain at least 
100 square resolution elements. Figure 26 shows the results of 
the sequential merging and the classifying process. In Figure 26 
there exist areas that were not classified therefore an additional 
pass was necessary to improve the classified areas. Figure 27 
shows small homogeneous areas clustered using the 36-square reso- 
lution element threshold. These areas are clusters 0, P, and Q. 
Figure 28 shows the areas 0, P, and Q, have been merged and classi- 
fied which improves the classified images. 


A. Data Problem Parameters Input Example 

$INPUT7 


NPASS = 2 
NCLUST = 0 

$END 

$INPUTA 


Perform two passes through the clustering 
merging and classifying process 5 

Initial cluster number. Additional passes 
can be performed on reruns, by inputing 
the total number of classes already detected. 


NSPS = 222 

NSCANS =120 
NCH = 12 
LT1 = 1 
LT9 = 9 

LT10 = 10 
LT11 = 11 
LT12 = 12 

LT13 = 13 

NSTART = 1 
NST0P = 140 


Number of resolution elements in entire 
scan or logical record 

Number of scans to process 

Total number of channels on input tape 

Logical tape unit for scratch tape 

Logical tape unit for storage of update 
information 

Logical tape unit for raw data input 

Logical tape unit for boundary map 

Logical tape unit for clustered homo- 
geneous areas 

Logical tape unit containing the final 
classified map 

Starting sample number 

Stopping sample number; Process only 
sample number 1-140 (NSTOP- NSTART + 1 255) 


NBTLG = 12 Bit length of word on input raw data tape 

MODE = 1 Input data FORTRAN formatted 

MODE = 2 non- FORTRAN formatted 
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Figure 25 Module Seven Output Example 
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Figure 28 Merged and Classified Homogeneous Areas 
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I TYPE - 0 
MSFC = 0 


14 = 1 


Input tape binary integer ITYPE ■ 1, 
floating point format 

Non-MSFC scanner format; MSFC scanner 
has housekeeping data requiring special 
handling. 

Set MSFC = 1 if MSFC scanner data 

Number of physical records per logical 
record 


NCRE = 1 

NSKIP = 0 
1NCX = 0 

1NCY = 8 

NSTX = 0 

NSTY - 0 


Data increment; NCRE = 1 use every data 
sample 

Initial records to skip before processing 

Increment in X direction on SC4020 plot 
frame for each sample 

Increment in Y direction on SC4020 plot 
frame for each sample 

Starting X coordinate on the SC4020 plot 
frame 

Starting Y coordinate on the SC4020 plot 
frame 


I XXX = 10 Homogeneous area threshold samples in X 

IYYY = 10 Homogeneous area threshold samples in Y 

direction 


© © 1234567890ABC ....(Alphanumeric characters for displaying 

classified homogeneous areas) 


B. Input Tapes 


Unit - A6 and B6 

A6 is any odd parity binary, 3- bit modulus, fixed 
point with word lengths <(_ 36 bits, or floating 
point binary. “ 

B6 is integer binary FORTRAN formatted. 

C. Output Tapes 

Units - B5 FORTRAN floating point binary 
(update information) 

A7 FORTRAN fixed point binary (Cluster data) 

B7 FORTRAN fixed point binary 


D. Intermediate Tapes 
Units - A3, Scratch tape only 
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Program Flow Chart 



set 

initial 

cluster 





Spectral Discrimination 








Spectral Discrimination 




Spectral Discrimination 















position 

tapes 



call 

FETCOR 



print 

^covariance 
\ matrix l 


print 

paper 


Spectral Discrimination 
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calculate 
eigen valves 
and vectors 



store class 
number in 
merg table 


assign 

population 


merge 

current 

cluster 
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i5 


Spectral Discrimination 








Spectral Discrimination 
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class 



1 


Spectral Discrimination 






119 








120 











/7/8 


// $DATA PROBLEM INPUT 


$IBFTC MOD 7B 


$ORIGIN B( MODULE 7B) 


/ _/ 

// $IBFTC MOD 7 A , 

^ORIGIN B( MODULE 7 A) > ‘ 


$IBFTC MOD 7 


$ORIGIN A{MODULE 7) 


= 3, 6 


/$IBFTC MODi 


/$ ORIGIN A(MODULEi) 


f $IBFTC MOD 2 


$ ORIGIN A{ MODULE 2) 


$IBFTC MODI 


/ $ ORIGIN A{MODULE 1) 


/$IBJOB FIOCS, MAP 


EXECUTE IB JOB 


$JOB 


Figure 29 IBM 7094 Job Deck Setup 
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IF EXCEEDS MAX: 


FAST tapes: a b c d 


□ STR QSTZ ODMPORETSY 


INPUT TAPES 



LINES OF OUTPUT (lOOO’S) MAXIMUW T,ME; 

DQ-i iA^-15 015-30 DOVER HOURS a3L MINI 
PROGRAMMER COMMENTS: [NUMBER OF CASES. 

CotY\ m e o + ( Cor* rrs e r\*\- 



OPERATOR COMMENTS: 


□ SEE ON-LINE 

O SEE TECHNIQUES 
n MAX EXCEEDED 
D RETURN TO SYS 

□ LINE MAX 


OPER injt:. 
over:- 


OUTPUT TAPES ONLY 




NO FILES NO FRAMES copies density 

1 fS ZIZZEl 


MSFC * Form 533 (Rev February 1966} 


Figure 30 IBM 7094 Instruction Form 533 
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APPENDIX 

Computer Program Listing 
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###*#*##**##*1 b M 7 O 0 & PROGRAM LISTING WITH JOB CARDS************* 


$ JOB NASA JAYROE BIN 53 *434300,00 , 12* I4MCE 

SEXECUTE I B JOB 

S 1 B JOB FI0CS,MAP 

$F ! LF -UNl T03- »NONF 

$F!LF -UNlTOA - »nonf 

$FILF -UN I T 07- * NONF 

SFILF -UNI T08-*N0NF 

SFILE -UNIT! 4 -* NONF 

SFILF -UNIT 1 5- * NONF 

SIBFTC FRFOP 

DIMENSION MODULE ( B ) 

COMMON DATAM?) 

namflist/init/modulf 

DO 1 1=1 ,8 

I MODULE (t 1=0 
1*0 

RFAD(5*INIT) 

WRITEJ6»INIT) 

498 CONTINUE 
T = T + 1 

IGO=MODULF( I 1 
IF (IGO .LF. 0) GO TO 999 
GO TO (10,11*12,13*14*15*161 *IGO 
10 CONTINUE 
CALL PROB 
GO TO 998 

II CONTINUE 
CALL GLEVEL 
GO TO 998 

1? CONTINUE 

CALL JONTPB 
GO TO 998 

13 CONTINUE 
CALL ISOMET 
GO TO 998 

14 CONTINUE 
CALL CNTLTN 
GO TO 998 

15 CONTINUE 
CALL BWNDR3 
GO TO 998 

16 CONTINUE 
CALL CLASFY 
GO TO 9©8 

999 CONTINUE 

CALL CLEAN 

STOP 

END 


SORIGIN ALPHA*SYSUT2*REW 
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/ 


SIBFTC MODI 

SURROUTINF PROP 
COMMON r>ATA< 1 7 ) 

COMMON/C1 / GLVLM 1 » 1 2 ) *NWH1CH< 1? J 
DIMENSION NT ABLE C 20 1 >1?) ,D4TAB( 201 ) 

NaMfL!ST/iNPUT1/Nch»NSPS*NScaNS»NSTaRT*NST0P,NbtLG*M00F» 
1ITYPF»MSFC*LTN,NSKiP, 
lNcRE*XMAX»XMiN,NOCHS*NWHICH 
RFADI S*INPUT1 ) 

WR I TF ( 6* I NPUT 1 ) 

C CALCULATE PROBILTTY DISTRIBUTION TABLE 

RFSOL=(XMaX~XMiN)/?O0.0 
SCALE=200.0/(XMAX-XMIN) 

C INITILIZE ARRAYS 

DO 6 I B= 1 *12 
DO 6 IBR=1»?01 
NTABLFt IPB»!R)=0 
6 CONTINUE 

DATAP ( 1 ) =XM I N 
DO 1 1=2*201 

DATAR ( I )=DATAB(I-1 l+RESOL 

1 CONTINUE 
NFLAG=0 

DO 4 I C= 1 *NSCANS 
NFLAG?=1 

DO 2 I=NSTART,NST0P,NCRE 

CALL GFT1 (DATA*NSP5*0,NCH*NSCAN0,LTN, |FRR»NFLAG2»NFLAG*NSTART* 
1NrTLG»MODF*NCRF» itypf*msfo 
GO TO (10*999*10,10,999) ♦ IFRR 
10 CONTINUE 

DO 3 I A = 1 * NOCHS 
NCHAN=NWHICH( IA) 

JJsDATA(NCHAN)*SCALF+l .0 

IF (JJ .LT. n GO TO 3 

TF ( JJ *GT. 20 1 ) GO TO 3 

NTARLEI JJ»NchaN)=NTABLF( JJ,NCHAN)+1 

3 CONTINUE 

2 CONTINUE 

4 CONTINUE 

999 CONTINUE 
REWIND LTN 
WRTTF(6*1000) 

1000 FOR MAT(1H1»SOX*26HPROR ABILITY DISTRIBUTION ) 

1001 FORMAT( 1H0,130h AMPLITUDE CH 1 CH 2 CH 3 CH 4 

1 CH 5 CH 6 CH 7 CH 8 CH 9 CH 10 CH 1 

71 CH 12 ) 

WRITE(6*1001 ) 

1002 FORMAT { 1H »FQ*1, 12110) 

DO 3 1=1 ,701 

WRlTF(6*l002) DATARt I) » (NTABLFt I , JJ) ♦ JJ=1 *NcH) 

5 CONTI NUF 

CALL PRODFNf NTABLF »SCALF»GLVL»XMAX ,XMTN> 

20 RFWIND LTN 
RETURN 
END 
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r> r\c\ 


*IBFTC PROAZ 

SUfiROUTlNP PROofN (NTapLf»SCALE»GLVL,XMAX»XMtN ) 
DIMENSION NTABLE ( 201 *12) *GLVL< 11*12) 

DO 3 |M*l ,12 

NTOTAL=0 

MAX=201 

DO 1 1=1, MAX 

NTOTALrNTOTAL+NTABLFU ,|M) 

1 CONTINUE 
NTFMPsO 

GLVLU*IM)=0.0 

nct=ntotal/io 

11=2 

DO 2 1=?, MAX 
I J=I-1 

NTFMP=NTFMP+CNTARLF(T J,IM)+NTABLF( I ,IM) )/2 
IF (NTFMP .LT. NCT) GO TO 2 
NTEMP=NTFMP-NCT 

GLVLU I, lM)=FLOATt I J ) /SCALF+XMI N-. S 
11 = 11+1 

2 CONTINUE 
GLVLtH ,1M)=XMAX 

3 CONTINUE 

DO 20 M = 1 ,11 

WRtTfI6»iO00) (GLVL(M,MM) »MM=1 *12) 

20 CONTINUE 

1000 FORMAT ( IX » 12E6. 1 ) 

RETURN 

END 


SIBFTC M0D2 

SUBROUTINE GLEVEL 

7094 DECK - MULT I SPECTRAL SCANNER PLOT 

DIMENSION I ALPHA ( 20 ) 

DIMENSION I AR AY ( 256 ) 

DIMENSION IWAVFM3) 

DIMENSION NUM(20) ,PeRCEN( 20) 

DIMENSION XLI M ( i 2 ) 

COMMON DAT A (1 7 ) 

C0MM0N/C1/ GL VL ( 1 1 * 1 2 ) »NWHICH(12) 

NAMeLTST/lNPuT2/NcH*NsPS*NsrANS*NsKTP,NSTART »NstOP , T TFRM,N, y CHAN , 

iiprt*iplt» iNcx,incy*nstx,nsty,nbtlg,mode,ncRe»ltn,iopt 
namelist /gRylvl/glvl 

I N=5 
I Nl =1 0 
l out =6 
NFLAGaO 
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n r\ r> -J A A 


nflgpt=o 

RFWIND INI 
CALL CAMRAV135) 

CALL RUTTV(l) 

RFAD( 5* INPUT 7 ) 

WRITFI6.INPUT?) 

00 303 1=1 *N 
XLTM(H=GLVL( I.ICH&N1 
303 CONTINUE 

RFAO (IN.GRYLVL) 

IF (10PT .FQ. 0) GO TO 302 
DO 301 I =1 » N 
XLlM(n=GLVLCI»ICHAN) 

301 CONTINUF 

302 CONTINUF 

********************* 

a-*#*#**#**#*******-#** 

RFAD(IN,1 02 ) 

RFAD (IN. 99) < I ALPHA < II . I = 1 .N ) .LLP 
RFAD( IN, 109 ) XLOW.XUPP* ( I WAVE ( I ) * 1 = 1.12) 

WR I TF ( I OUT » 11 0 ) 

1003 FORMAT (1X»12A6»16»5X»2F8«4) 

300 CONTINUF 

00 777 1=1 » N 
NIJM(I)=0 
77 CONTINUF 
KK=N-1 

WRlTFdOUT.106) 

WRITE( I OUT . 1 10 ) 

WRITFUOUT.102) 

WRTTFf IOUT.106) 

WR 1 TF I TOUT *1031 TCHAN 
00 3 1=1 «KK 
I F ( 1 • F0« 1 1 GO TO 9 

WRlTF (6.105) I ALPHA( I) »XL!M{ T ) *XLIM{ I+H 
GO TO 2 

9 WR|TF( I0UT.112) ( IWAVF(K) *K=1 » 12 > * I ALPHA ( I ) .XLlM ( I ) ,XL IM (1+1 ) 

GO TO 2 
2 CONTINUF 

WRTTFUOUT.106) 

CALL LARFLl ( NST aRT , NSTOP ♦ N rRF ) 

1 NUM COUNTS SCAN LINES 


INUM=0 

IDIFsNSTOP-NSTART+1 
IF (NSKIP ,F0* 0) GO TO 806 
00 805 1=1 .NSKIP 
NFL AG3=t 

CALL GFT1 (OATA,NSPS,0,NCH.NSCANO,LTN,IFRR.NFLAG2.NFLAG.NSTART, 
1NRTLG.M0DF.NCRF. I TYPE. MSFC) 

805 CONTINUE 

806 CONTINUF 
NFLAG1 =0 

8 CONTINUF 
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IINaO 

NFLAG?=1 

205 00 90 U=NSTaRT*NST0P 

!IN«!IN+1 

CALL GET1 (DATAtNSPS* 0 ♦ NCH .NSCANO ♦ I Nl ♦ I ERR *NfLAG2 »NfLAG »N$TaRT , 

1N0TLG*«0DP»NcRF» ITYPF*MSFC) 

GO TO(800»R01 .R00,R00*801 ) ,IFRR 
801 WR ! TF l TOUT ♦ 1 00 ) NSCANO 
GO TO 900 
800 CONTINUF 

DO 601 Is?*N 

IF (DATA(ICHAN) «GT. XLIM(I)) GO TO 601 
NuM( I~ 1 ) = NUM( 1 - 1 >+1 
IARAYUIN) = IALPHAU-1 ) 

GO TO 90 
601 CONTI NUF 
90 CONTINUF 

numtot=inum*idif 

totnum^numtot 

TOTPFR=0.0 
DO 880 1=1 » 11 
XXXNUMrNUMU ) 

PFRCPNt n=(XXXNUM/TOTNUMl * 100.0 
TOTPFR=TOTPFR+PFRCFN ( I ) 

888 CONTINUF 
96 CONTINUF 

INUM=INUM+1 

IF (IPRT .FQ. 0) GO TO 889 

WRTTFt IOUT»10A1 NSCANO, UARAY(K) ,K=l , IDIF) »LLP 

889 CONTINUF 

IF UPLT .FQ* 0> GO TO 890 

CALL PLT8F1 UARAY, 10IF»NBLK,INCX* INCY»NSTX,NSTY,NCRF. 
1NfLGPT*NfLAG1 ) 

890 CONTINUF 

IF (NSCANO .GF. NSCANS) GO TO 218 
GO TO 5 
218 CONTINUE 
NFL AG1 =0 
NFLGPT = 0 
NSCANO=0 
N5TART=NST0P+I 
NSTOPsNSTOP+IDIF 
IF (NSTOP .GF. NSPS) NSTOP=NSPS 
ITFRM=ITFRM-1 
REWIND INI 
WRITE CIOUT.llO) 

WRITE ( I OUT *114) 

DO 666 1*1 *11 

WRlTP ( lOUT *115) I ,NUM( n *PFRCEN( I ) ,TALP hA( I } 

666 CONTINUF 

WRTTF ( I OUT » 1 16 ) NUMTOT.TOTPER 
1F(ITFRH,GT*0) go TO 300 

98 WRITF(IOUT*HO) 

99 FORMAT ( 25A1 ) 

100 FORMAT (161 4) 
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101 FORMAT ( 12F6.1 ) 

10? FORMAT (flOH 

1 > 

103 FORMAT (2X*8HCHANNEL I 2 *36X * 22HSYMBOL INTERVAL /) 

104 FORMAT (5X* 14* IX » 1H* * 1 20A1 *Al ) 

105 FORMAT< 50X *A1 *8X*F6.1 *3H - *F6.1) 

10ft FORMAT ( 1 X / ) 

105 FORMAT(lOX*115Al ) 

109 FORMaT(2F10.0*12A1 ) 

no formatiihi/i 

111 FORMATC2X.83HCROPS C5-17 *5-22* 5-5 , 5-8* ALTITUDE 2000 FT. 

1 DATE RFCOROFD 6/30/66 //) 

112 FORMAT C2X»23hWaVF LFNgTH (MICRONS) * 12 Al » 13X * Al »8X » F6 , 1 * 

1 3H - *F6 • 1 ) 

113 F0RMAT(nx»112Al ) 

114 f0RMaT(20X*5hCLaSS,7X,13HSAMPLF NuMrfR * 6X *7hPfRceNT * 10X » 

1 12HALPHA SYMrOL/) 

115 FORMAT (21X*i2»10X*t6»12X»f7.2»15X»a1) 

116 FORMAT l 20X * 5HT0TAL * 8X » 1 6 » 1 2X * F7 . 2 ) 

117 FORMAT <3I4»12F7.3) 

999 CONTINUE 

REWIND INI 

return 

END 


SIBFTC PL0TT1 

SURROUTINF LA8FLHNSTART*NST0P*INCRE ) 

DIMENSION I0UT( 120) ,JOUT( 120) *<OUT( 120) ♦ LOUT ( 120) 
NDIF=(NST0P-NSTART+1 J/INCRE 
11*0 

DO 1 I=NSTART.NSTOP,lNcRE 
!!=!I+! 

lOUTdl ) = I/1000 
JOUTdl ) = I/1 00^1/ 1000*10 
KOUT (II )=I /10- 1/100*10 
LOUTdl ) = I-I/10*10 
IF (LOUTdl) .LE. 0 ) LOUT ( 1 1 ) =0 
1 CONTINUE 

WRtTf (6*10) dOUTd ),I = 1 *NdIF) 

WRlTF (6*10). (JOUTd ) *1 = 1 *NDTF) 

WRlTF (6*10) ( KOUT (I)*I=1*NDTF) 

WRlTF (6*10) (LOUT ( i ) » 1=1 »ND!F ) 

io format (nx*i20in 

RETURN 

END 


5IRFTC PLLTT1 

SURROUTINF PLTBF1 ( I SYM ,NN , NsQR ♦ I NcX » I NcY *NSTX ♦ NSTY ,N cRf» 
1NFLAG»NFLAG1 ) 

DIMENSION ISYM(1 ) *NRUFER(50) 

IF (NFLAG .NF, 0) GO TO 10 
NSOR*10?4/IABS(INCY+INCX) 

10 CONTINUE 
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IF (NFLAG1 »NE, 0) GO TO 20 
CALL pRAMFV(O) 

NCOUNT=0 

INcRXsNSTX 

inc r y=nsty 

NFLAG1=1 

IF (NFLAG .FQ. 0) GO TO 11 

CALL APRNTV ( INCX.INCY*NN,N(3UFFR* INcRX»INcRY) 
INCRYalNcRY+lNCX 
I NCRX = 1 NCRX+t ABS ( INrY ) 

11 NFLAG=1 
20 CONTINUE 

NCOUNTsNCOUNT+i 
DO 1 1=1,50 
N8UFFR (11=0 

1 CONTINUE 
I A = 0 

DO 2 1 = 1 , NN 
I B= I A+6 
I C= ! B /6 
ID=IB-IC*6 
IE=ID*6 
I A=I A+l 

CALL FLDtNRUFER( IC ) »IF»6» OtlSYM(l)) 

2 CONTINUE 

CALL aPRNTV ( I N cX , I NrY ,NN ,NpUFE R » T Nc R X , I N^RY ) 
INCRY=INCRY+INCX 
I NC r X = I NC R X+ I AR5( INCY ) 

IF CNCOUNT .GE, NSOR) NFLAG1=0 

RETURN 

END 


SIBFTC GFTDA1 

SUrROUTiNe GET1(DATA,NSPS,NsKiP,NcPw,NSCAN0, iRw, IE rr ,NfLA62, 
INFLAG* 

1NSTART,NrTLg»M0dE,NcRE» ITYPe»MSEC ) 

DIMENSION DATA ( 1 ) »NoAT { 890 j 
DATA NFLAG/O/ 

IF (NFLAG .NE. 0 ) GO TO 12 

ntemp=ncpw*<ncre-i ) 

NSCANOsO 

NL*36-N8TLG 

NBLNG=36/NBTLG 

NFLAG=1 

NSTMP=(NSPS+NSTART-l ) *NCPW+MSFC*19 
NSCANS=NSTMP 

ntmp*nstmp/nblng 

nstmp=nstmp/nblng*nrlng 

if (NSTMP .LT, NSCANS) NTMP=NTMP+l 
NSCANS=NTMP 

IF (NSKIP • EQ • 0) GO TO 402 
DO 401 I=1»NSKIP 
NSCAN0=NSCAN0+1 
CALL SXRBIN < IRW,1 »RD } 
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401 CONTI NUF 
40? CONTINUE 
1? CONTINUE 

IF (NFLAG2 .FQ. 0) GO TO 10 

13 CONTINUE 

Ma(NSTART-l )*NCPW+MSFC*19 
IF (MODE .EQ. 1) GO TO 50 

CALL RFDTPR ( IRW.MODF. IFRR .NW.NSCANS.NDAT } 
GO TO 51 

30 CALL REOTPC ( iRW .MODE » IFRR .NW.NS caNS.NdAT ) 

51 CONTINUE 

1000 FORMAT < 1 X • 1 6 > 

NNW=NW*6 
NFLAG2=0 
NSCANOsNSCANO+1 
10 CONTINUE 

DO 14 NN = 1 » NCPW 
IRrM+NBLNG 

id*ib/nblng 

if=ib-nblng*id 

ibit=nbtlg*if 

MaM+1 

IF CITYPP .EG. 1) GO TO 15 
NDAT A = 0 

CALL FLD (NDATA.NL .NBTLG* IB I T»NDAT( ID)) 
DATA(NN)=NDATA 
GO TO 14 
15 CONTINUE 
DATANsO. 0 

call FLd( DATaN.NL .NbTLG* IB IT »NDAT( ID)) 

DATA (NN ) =DATAN 

14 CONTINUE 
M=M+NTEMP 
RETURN 
END 


SORIGIN ALPHA. SYSUT2.REW 

SIBFTC M0D3 

subroutine JONTPB 

DIMENSION NKNT<4000) .NP(4000.2) »NWKCD( 750).NWKcL( 750) 
DIMENSION DATAM?) 

C DIMENSION I B I N ( 750) 

DIMENSION ALPNUM ( 47 ) .ALPHA (120) *CORDX (3) 

DIMENSION IMX(6) »IMY(6> 

NAMELIST /INPUT3/ NrH .NSPS .NSCANS »NS< I P . NST aRT .NSTOP ,LtN .NbtLg. 
1MOdE.ITYPe.MSFC.NCRF.NOJP. I MX. I my. SCALE. B I AS 
DATA ASTRIK/1H*/ 

DATA XMARK/IHX/ 

DATA RLANK/6H / 

1060 FORMAT (48X.25HDATA SWITCH HAS OCCURRED ) 

1061 FORMAT ( 49X . 30H JO I NT PROBABILITY DISTRIBUTION ) 
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1062 FORMaT(IHI) 

1063 formatiaax.iih x-axts is » 1 6 »6x 1 1 ih y-axis is ,i6) 

1040 FORMATdH *67H MAXIMUM PROBABILITY OF L'NcOMMONaI I TY EXCEEDED” CONt 
1INUE EXCUT I ON ,16) 

1064 FORMAT ( IX ♦26HSYMB0L N/SYMBOL ) 

1065 fORMaT(11X»121 (1h*1*/*UX*1h*»55X,5hPaRT ,I1,4h0f * 1 1 . 53X » 1H* . / 

1 • 1 1 X * 1 2 1. ( 1 H* ) ) 

IRT=1 1 
NFLG3=0 
LIM|T=4000 
READ( 5* INPUTS ) 

WR|TE(6*INPUT?) 

RFAD ( 5 » 1000 ) (ALPNUMf 1 ) *1 = 1 .47) 

1000 FORMAT ( IX » 47A1 ) ‘ 

C RFWIND IRT 

DO 80 NT IMF=1 *NOJP 

NFL6=0 

NH»1 

NV=2 

IF (NSKIP ,E0. 0) GO TO 9 

DO 8 I=1,NSKIP 

CALL SKRBINILTN.1 ,RD) 

8 CONTINUE 

9 CONTINUE 
MIX*IMX(NTIME) 

MIY=IMY(NTIME> 

DO 10 1=1.4000 

10 NKNT ( I ) =1 
WRITE ( 6 » 1 062 ) 

WRITFI6.106I ) 

WRTTF(6*l063) MIX.MTY 
NCLS=1 

DO IB MRP= 1 .NSCANS 
NFLAG2=1 

DO 15 MPR=1,NSPS 
IF (NFLG3 .GT. 0) GO TO 16 

CALL GFT2 tDATA.NSPS.0.NCH»N5CAN0 ,LTN, ieRR»NfLaG2.NSTaRT» 

inbtlg*modf.ncRf» itype.msfci 

Nc=DATA(MlX)*SCALE+oi as 
NR=DATA(MIY)*SCALE+RIAS 
tF (NC .LT. 1) NC=1 
IF (NC .GT. 255) NC=255 
IF (NR .LT. 1) NR=l 
IF (NR ,GT. 255) NR=:255 
NP ( 1 ,NH)=NC 
NP( 1 *NV ) =NR 
NFLGBsl 
GO TO 15 
16 CONTINUE 

CALL GET2 (DATA.NSPS.O.NCH.NSCANO.LTN, ifRR,NfLAG2.NsTaRT. 
lNBTLG»MODEtNcRF. I TYPE .MSEC) 

GO TO (996.999.998.996.999 ) , IERR 
999 WRITF(6.1 070) NSCANO 

1070 FORMAT( 33HFNO-OF-FILF ON INPUT AT SCAN NO. .16) 

GO TO 995 
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998 WRITF(6.1 071 ) 

1071 FORMAT<27HIRRECOVERABLE PARITY ERROR > 

60 TO 995 
996 CONTINUE 

NCaOATAIMlX )*SCALF+RI AS 
NR=DATA(MIY)*SCALF+RI AS 
IF (NC .LT. 1 ) Nr=l 

IF (NC .GT • 255) Nr=255 

IF (NR .LT. 1) NR= 1 

IF (NR .6T. 255) NR=255 

DO 11 1=1 *NCLS 

IF INC .NE. NP(I«NH)) GO TO H 

IF (NR .NE. NP(I.NV)) GO TO U 

NKNT ( I ) *NKNT ( I )+l 
GO TO 15 
11 CONTINUE 

IF (NCLS .GE. LIMIT) GO TO 15 
NCLS=NCLS+1 
NP(NCLS.NH)=NC 
NP(NCLS*NV)=NR 
15 CONTINUE 
13 CONTINUE 
17 CONTINUE 
1020 FORMAT (IX, 316) 

C FIND MINIMUM AND MAXIMUM IN ROW AND COLUMN DATA 

Ml NC=NP ( 1 » NH ) 

MAXC=NP( 1 »NH) 

Ml NR=NP ( l * NV ) 
maxr*np< 1 »NV ) 

DO ?5 1=1 » NCLS 

IF (NP< I *NH) .GT. MAXC) MAXC=NP (I »NH ) 

IF ( NP (I * NV ) .GT. MAXR) MaXR=NP( I *NV) 

IF ( NP ( I * NH ) .LT, MlNC) MINC=NP ( I »NH) 

IF ( NP ( I * NV ) .LT. MINR) M I NR = NP ( I »NV ) 

25 CONTINUE 
C TEST FOR SPREAD 

ndifv=maxr-minr+i 

ndtfh=maxc-minc+i 

IF (NDIFV .GE. NDIFH1 GO TO 19 
WR I TP { 6 ♦ 1 060) 

WR I Te ( 6 * 1 063 ) Ml Y*MlX 

NH= 2 

NV=1 

NTEMP=MAXC 

MAXCsMAXR 

MAXRsNTEMP 

NTFMPsMINC 

MI NC*MI NR 

MINR=NTEMP 

NTEMPsNDIFH 

ndifh=ndifv 
ndifv=ntemp 
19 CONTINUE 

WRITF(6*1 064) 

C ORDER ROW DATA IN 0 F CENDI NG ORDER 
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n r\ 


C BUBBLE-UP 

NSORTrNCLS 

21 CONTINUE 

DO 20 1=2 »N50RT 

IF (NP(I-1*NV) ,GE, NP(I,NV>) GO TO 20 
NTFMPl=NP{I-i.NV) 

NT6MP2=NP ( I-i *NH ) 

NTFMP^=NKNT( I-l > 

NP{ I-l »NV)=NP( I *NV ) 

NP ( I-l »NH)=NP( I »NH) 

NKNT( I-l ) =NKNT ( I ) 

NP C I *NV ) «NTEMP 1 
NP ( I ,NH)=NTFMP? 

NKNT( IJ=NTEMP3 ‘ 

20 CONTI NUF 

NSORT=NSORT-l 

IF (NSORT .GT* 1) GO TO 21 

C calculate tablf 

MAXKNT=0 
DO 22 1=1 >NCLS 

22 IF (NKNT(I) # GT. MaXKNT) MaXKNT=NKNT < I ) 
NFACT=MAXKNT/45 

XXX*FLOaT l MaXKNT ) /A6.0 
NFAC=0 

IF (NFACT ,LT* 1) NFACT=1 
WRlTF(6»l050) BLANK, NpaC 
NFAC=NFAC+NFACT 
DO 23 1=1 >46 

WRlTF(6>1 050) ALPNUMU),NFAC 
1050 FORMAT (3X,A6>6X,I6) 

NFAC=NFAC+NFACT 

23 CONTINUE , 

GO TO 18 

14 CONTINUE 

WR I TE 1 6 > 1 040 ) LIMIT 
GO TO 17 
18 CONTINUE 

PRINT DISTRIBUTION ON PAGE 

SET UP LOOP 

ITFRM=1 

IF (NDIFH .EQ. 120) GO TO 40 
ITFRM=NDIFH/120+1 
40 CONTINUE 
NST * 1 

NSO=NST+l 19 

IF (NSO .GT. NDIFH) NSO=NDIFH 
DO 50 LOOP = l » I TERM 
WRITE(6>1065 ) LOOP , I T FRM 
CORDXI1 )=FLOAT (NST+MINC-2) 
CORDX(2)=COROX( 1 ) +60 • 0 
CORDX(3)=CORDX( 1 ) + 11 0.0 
1021 FORMAT ( 1 H 1 ) 

I B = 1 

DO 60 II=1»NDIFV 
DO 66 1=1 >120 
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ALPHAd 1=BLANK 

66 CONTINUE 

C IF (NFLG «EO • 1) GO TO 69 

C DO 68 1=1 »NOIFH 

C I BI N ( I ) =0 

C68 CONTINUE 
C69 CONTINUE 

IF (NP(IR.NV) .NE. MAXR-TI+l) GO TO 68 
I A=0 

00 26 UIBtNCLS 

1 F= I 

IF CNP( I *NV ) ,NE. NP(IB»NV)) GO TO 31 
I A* I A+l 

NWKCDl I A ) =NP ( I »N'H) 

NWKCL(IA)=NKNT(I } 

26 CONTINUE 
31 CONTINUE 

I B* 1 E 

C ORDER X-AXIS 

NSORT= I A 
NUMCLX=IA 

IF ( IA .LE. 1 > GO TO 30 

28 CONTINUE 

DO 27 IC=?»NSORT 

IF CNWKCDt IC-1 ) ,LE* NWKCDUC)) GO TO 27 
NTEMP1=NWKCD( IC-1 > 

NTEMP2=NWKCL( IC-1 ) 

NWKCD ( IC“1 >=NW<CD< TO 
NWKCLUC-l )=NWKCL( ID 
NWKCD(IO=NTFMPi 
NWKCL(IC)=NTFMP? 

27 CONTINUE 
NSORT=NSORT-1 

IF ( NSORT .GE. 2) GO TO 28 
1030 F0RMATI1X .216) 

30 CONTINUE 

C IF (NFLG .EQ. 1 ) GO TO 91 

C IM*1 

C DO 90 1 = 1 .NDIFH 

C IF (NWKCD(IM) .NF. T> GO TO 90 

C I BI N ( I ) =NWKCL < I M ) 

C I Ms I M+l 

C90 CONTINUE 

C91 CONTINUE 

NMl=NST+MINOt 

DO 67 I D=1 .NUMCLX 

IF (NWKCD(ID) ,LT. NM1 ) GO TO 67 

IC=T0 

GO TO 29 

67 CONTINUE 

29 CONTINUE 
I A = 0 

DO 64 I =NST »N50 
I A= I A+l 

IF (NWKCD(IC) *NF • MINC+I-I) GO TO 64 
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XXsFLOaT (NWKcU IO ) /XXX 
1CaR*XX+(1. 001-1. O/XXX) 

ALPHA<!A)*ALPNUM( ICAR) 

ic=ic+i 

64 CONTINUE 

65 CONTINUE 

IP (NFLG • EO. 1) CO TO 70 
C WRTTFURT) ( IBINf I ) *1=1 *NDIFH) 

70 CONTINUE 

CORDY=FLOAT (MAXR-I 1+1 ) 

YMARGsXMARK 

IF (NOIFV-II .NE. (NOIFV-1-!I/10*10)| YMaRG=ASTRiK 
WRlTFt6»l008) cORhY*YMaRg* ( ALPHA ( I) *1=1*120) 

1008 FORMAT <1X*F8.1*2X*A1 *120 A 1 ) 

60 CONTINUE 
NST=NSO+l 
NSO=NST+l 1 9 

IF CNSO ,GT. NDIFH) NSO=NDT FH 
WRITEC6»1010) 

WRiTf(6*10id CORDXII ) »C0RDX{2) *CORdX(3) 

1011 FORMAT {6X*F10.4,50X*F10.4*40X*F10.4) 

NFLG*1 

WRlTF(6*l06?) 

50 CONTINUE 

REWIND LTN 
C END FILE IRT 

80 CONTINUE 
C REWIND IRT 

995 CONTINUE 

10)0 FORMAT M1X*1?(10hX*********)) 

RETURN 

END 


*IBFTC GETOA2 

SUBROUTINE GET?{DATA»NSPS*NSKIP,NcPW*NSCANO, iRw* IERR,NfLAG 2* 
lNSTART.NRTLG»MODE*NrRF*lTYPF»MSFC > 

DIMENSION DATAU ) »NDAT(890) 
data NFLAG/O/ 

IF (NFLAG .NF. 0 ) GO TO 12 
NTEMP=NCPW*(NCRE-1 ) 

NSCANO=0 

NL=36-NBTLG 

NBLNG=36/NBTLG 

NFLAGal 

NSTMP=(NSP$+NSTART-1)*NCPW+MSFC*19 

NSCANS=NSTMP 

NTMP=NSTMP/NBLNG 

nstmp=nstmp/nblng*nblng 

if (NSTMP .LT, NSCANS) NTMP=NTMP+1 
NSCANS=NTMP 

IF (NSKIP .EQ. 0) GO TO 402 
DO 401 1=1 *NSKIP 
NSCANO=NSCANO+i 
CALL SKRRIN (IRW*1,RD ) 
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401 CONTINUF 

402 CONTJNUF 
12 CONTINUE 

IF (NFLAG? .EQ. 0) 60 TO 10 
16 CONTINUF 

M*{NSTaRT- 1 )*NcPW+MSFC*19 
IF ( MODF ♦ EQ. 1 ) GO TO 50 

CALL RFDTPR ( lRW»MODF*If r RR*NW*NSCANS»NOAT) 
GO TO 51 

50 CALL REDTPC < lRW*MOOE» lERR *NW»NSCAN$,NDAT ) 

51 CONTINUE 

1000 FORMAT ( 1 X * 1 6 ) 

NNW=NW*6 
NFLAG2=0 
NSCAN0=NSCAN0+1 
10 CONTINUE 

DO 14 NN=1»NCPW 
IB=M+NBLNG 

id=ib/nblng 

IF»IB-NBLNG*ID 

IBIT=NBTLG*IF 

MrM+1 

IF ( ITYPF .EQ. 1 > GO TO 15 
NDATA=0 

CALL FLD(NDATA»NL»NBTLG*IBIT»NDAT( ID) ) 
DATA(NN)=NDATA 
GO TO 14 
15 CONTINUE 
DATAN=0. 0 

CALL FLD(DATAN,NL»NBTLG» IBIT»NDATUD) ) 
DATA(NN)=DATAN 
14 CONTINUE 
MsM+NTEMP 
RETURN 
END 


SORIGTN ALPHAiSYSUT2*REW 

SIBFTC M0D4 

subRoutjnf ISOMFT 

DIMENSION NTA<256) .NTR(256) ♦ DATA I 256 ) * T I MF ( 256 ) *DAT ( 12 ) 

DIMENSION NTCI256) »NTD(256) 

namelist/ i npuT4/ nch*nsPs.nskip,nbtlg,mode*iRw»nchan,nsncRe,npcRE* 

1 1 TYPE »MSFC ♦ 

2NPTSL*NPTSU, 

1 MaXScN*XMi n ♦xMax*ymin»ymax*nblszx»nblszy*nsect* nsmov 
?*NDIREC 
READ t 5 » I NPUT 4 ) 

WR t TE ( 6 * I NPUT 4 ) 

NUP*5 1 2/NBL5ZY 

nbldiR=nblszx*ndtRec 

ntimfs=maxscn/nup 

CALL CAMRAV(O) 

nback=i 
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IF ( NSK I P *F0* 0) GO TO 33 

DO 32 1=1 »NSK I P 

NFLAG2M 

CALL GFT 3 ( DAT »NSPS. 0 .NCH »N$CANO, i Rw. I fRR ,N pLAG2 .NPTSL .NBTLG. 
12 *NPcRF» ITYPE»MSFC) 

32 CONTINUE 

33 CONTINUE 

DO lfl MMRR=1 *NSECT 

NPTSUU=NPTSU 

XMINsFLOAT(NPTSL) 

XMAXsFLOAT(NPTSUU) 

NDIFF=NPTSUU-NPTSL+t 

IF (NSMOV .Np* 0) NPT$UM=N0IFF -NSMOV+1 
nflag=o 

DO 16 MPR=1 »NT IMFS 
IF (NDIRFC) 23*24.24 

24 IXR=512-( 512/NrLSZX-NUP)*NrLSZX 
IXL=0 

GO TO 25 

23 IXL=512“( 512/NBLSZX~NUP>*NBLSZX 

I XR = 0 

25 CONTINUE 
IYR=512 

I YT = 0 

CALL FRAMEV(O) 

DO 10 MRP=l.NUP 
DO 17 MMRP=1 .NSNCRE 
NFLAG2=1 

CALL GFT 3 ( DAT .NSPS.O.NCH.NSCANO, I Rw. IERR.NFLAG2 .NPtSL .NBTLG. 

i mode *npcRe*i type. m$fo 

IF (NFL AG *FQ* 0) GO TO 19 
17 CONTINUE 

19 CONTINUE 
NFLAG=1 

DATA(1)=DAT(NCHAN) 

DO 1 1=2. ND IFF 

CALL GET3(DAT .NsPS, 0 »NcH.NsCANO,lRw.IERR.NFLAG 2 .M p TSL,NBTLG. 
imode.npcRe* ITYPf.MSFC) 

DATA ( I 1=DAT(NCHAN) 

1 CONTINUE 

IF (NSMOV .EQ. 0) GO TO 22 
DO 20 1=1 .NPTSUM 
I A — I 

IB«I a+nsmov-i 
temp=o.o 

DO 21 1 1 = I A+ 1 B 
TFMP=TEMP+DATA(1 I ) 

21 CONTINUF 
TFMP=TEMP/FLOAT ( NSMOV } 

DATA ( I ) =TFMP 

20 CONTINUF 

22 CONTINUF 

CALL SETMIV ( I XL » I XR . I YB . I YT ) 

CALL XSCALV (XMlN.XMAX.IXL.tXR) 

CALL YSCALV (YMlN.YMAX.tYR.lYT ) 
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00 11 1=1 »NDIFF 

T|ME( 1 ) =FLOaT ( NPTSL- l+I ) 

NTACI ) *NYV ( DAT A ( I ) ) 

NTR( n=NXV(TTMF< I ) ) 

11 CONTINUE 

DO 12 I =2 *NOI FF 

CALL LINFV (NTB( 1-1 ) »NTA(I-1 ) »NTB ( I )*NTA<I )) 

1? CONTINUE 

IF (MRP .LE. 1 ) GO TO 15 
00 15 1=1 *NDI FF 

CALL LINFV (NTDdl.NTCd) »NTP(n*NTA(IH 

13 CONTINUE 

15 CONTINUE 

DO 1A 1=1 *Noi FF 
NTOC I)=NTB( I J 
NTC ( I ) =NTA ( I ) 

14 CONTINUE 
IXL=IXL+NBLDIR 
IXR=IXR-NBLDIR 
I YB= I YB-NBLSZ Y 
IYT=IYT+NBLS2Y 

10 CONTINUE 

CALL BSRFCDt I RW » NBArK ,NOF ) 

16 CONTINUE 
rewind irw 

IF (NSKIP .EQ. 0) GO TO 30 
DO 31 1=1 *NSK I P 
NFL AG2=1 

CALL GET 3 (DAT .NSPS , 0 ♦ Nch »NSCANO , I Rw » I fRR * NfL AG2 »NPtSL ,NbTLG. 
1M0dE*NPcRE»ITYPF,M5FC) 

31 CONTINUE 
30 CONTINUE 

nptsu=nptsuu 

ntemp=nptsl 

nptsl=nptsu+i 

NPTSU=NPTSU+(NPTSU-NTEMP )+i 
18 CONTINUE 
RETURN 
END 


SIBFTC GET0A3 

SUBROUTINE GET3(DATA»N$PS,NSKIP,NCPW.NSCAN0, IRW,IERR,NFLAG2. 

1 nst art ♦ nbtlg » mode *ncRe» itypf»msfc> 

DIMENSION DATA(1 ) *NDAT(890) 

DATA NFLAG/O/ 

IF (NFLAG .NF. 0 ) pO TO 12 

NTEMP=(NCRE-1 )*NCPW 

NSCANOsO 

NL=36-N8TLG 

NBLNG=36/NBTLG 

NFLAG-1 

NSTMP=(NSPS+NSTART-1 )*NCPW+MSFC*19 
NSCANS=NSTMP 

NTMPsNSTMP/NBLNG 
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nstmp*nstmp/nblng*nrlng 

if < NSTMP .LT , NSCANS) NTMP=NTMP+1 
NSCANS»NTMP 

IF (NSKIP .EQ. 0) GO TO 40? 

DO 401 I=1»NSKIP 

NSCAN0=NSCAN0+1 

CALL SKRRIN ( I RW * 1 » RD ) 

401 CONTINUE 

402 CONTINUE 
1? CONTINUE 

IF (NFLAG2 .EQ. 0) GO TO 10 

13 CONTINUE 

Ms(NSTART“l ) *N cPW+M5FC*1 9 
IF (MODF .EQ* 1) GO TO 50 

CALL REDTPR (IRW*M0DE»IERR*NW.NSCANS.NDAT) 
GO TO 51 

50 CALL REDTPC ( IRW.MOdE » I ERR ^^SCANS » NDAT ) 

51 CONTINUE 

1000 FORMAT ( 1 X ♦ 1 6 ) 

NNW=NW*6 
NFLAG?=0 
NSCAN0=NSCAN0+1 
10 CONTINUE 

DO 14 NN=1.NCPW 

ib=m+nblng 

io=ib/nblng 

if=ib-nblng*id 

ibit=nbtlg*if 

M=M+1 

IF ( I TYPE .EQ. 1) GO TO 15 
NDATAsO 

CALL FLDl NDATA.NL. NBTLG* IB IT *ND AT (ID) ) 
DATA(NN)=NDATA 
GO TO 14 
15 CONTINUE 
DATAN=0.0 

CALL FLD(DATaN,NL*NrTLG»IBIT .NDAT ( ID) ) 
DATA(NN)=DATAN 

14 CONTINUE 
M*M+NTEMP 
RETURN 
END 


SORIGIN ALPHA.SYSUT2.REW 

SIBFTC MODS 

subroutine cntlin 
DIMENSION DAT(1?) 

COMMON /Lab /Ni3»Nf3.NT3.Z{2560)*KON(3400)»Ni2.Nf2.NT2.U(3400)» 
1 VC 3400) 
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namelist /inputs/ nch,nsPs.irw»nchan,n$ncre ,npcre*nptsl,nptsu* 

2NBTLG*M0DE* 

2 1 TYPE »MSFC » 

3NSKIP* 

lMAXScN*NSFCT» MSZX. MSZY* RLK,FHINc»ZM!N*ZMaX*LaR 
RFAO (5*1 NPUT5 ) 

WRITE (6*1 NPUT 5 ) 

CALL OPEN ( IRW* 1 ) 

lXR=l024.0-BLK+.5 

IXL*0 

IYR = 0 

I YT = 1 36 

NCRE=NPCRE 

NBLLK=l024*0/BLK+.5 

XNBLKs blk 

FNPTSL=NPTSL 

FNPT SU=NPT SU 

NUPsMAXSCN/ ( (MSZY-l )*NSNCRE) 

NSECT=NSPS/ MSZX 
NDIF=(NPTSU-NPTSL+l )/NPcRE 

NFLAG1=0 
FDIF=NDIF 
FRM= MSZY 
CALL CAMR AV ( 35) 

NWDS=NCH*NSPS 
CALL BUTT V < 1 ) 

CALL SFTM]V( IXL*IXR*IYR*IYT) 

IF (NSKIP • EQ« 0) GO TO 25 
DO 26 1=1 *NSK I P 
NFLAG2=1 

call get5(dat*nsps *o*nch*nscano,irw» IERR*NfLA62»NPTSL* 
1NbTLG*M0dE*NPcRE* ITYPe*MSFC) 

26 CONTINUE 
25 CONTINUE 
C READ(IRW) 

1 1 1 = 0 

CALL FRAMFV(O) 

DO 18 MMRR=1*NSECT 
NPLTMPsNPTSL-f 1 

nptsuu=nptsu 

NPTSUsNPTSU-NPcRE 
NFLAG=0 
I MM=0 

DO 16 MRPP=l,NUP 
111 = 111+1 

IF (III .LE. NBLLKJGO TO 21 
CALL FRAMEV(O) 

111 = 1 

lXR=l024.0-BLK+*5 
I XL*0 

21 CONTINUE 

CALL YSCALVU.O ♦ EDI F »IYB»IYT> 

CALL XSCALV<1«0»fRM*IXL*IXR) 

DO 10 MRP s 1 » MSZY 
IF (NFLAG1 *NE« 1) GO TO 20 
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n n 


DO 22 I1M=1,NDIF 
I MM = IMM+1 
IMI*IMI+1 
Z( TMM>=Z( IMJ ) 

22 CONTINUE 
NFLA61=0 
60 TO 10 
20 CONTINUE 

DO 17 MMRP=1 *NSNCRE 
NFLAG?=1 

1001 F0RMaT(1X»16F4»1 1 

CALL GET5<DAT»NSPS »0 »NcH »NSCANO » I R w » IERR *N fLaG2 »NPt$L , 
1NbTLG*M0DF»NPcRF* ITYPF»MSFC) 

IF (NFLAG .EQ. 0) GO TO 19 

17 CONTINUE 
19 CONTINUE 

NFLAGM 

I MM = I MM+1 

Z ( I MM ) =OAT ( NCHAN ) 

DO 1 I=NPLTMP,NPTSUU*NPcRE 
IMMsIMM+1 

CALL GET5(DAT*NSPS * 0 *Nch ,NSCANO* \ Rw , I ERR *NpLAG2 .NPjsL , 
1NbTLG*M0DF»NPcRE»ITYPF*MSFC) 

Z( I MM ) =DAT ( NCHAN ) 

1 CONTINUE 
10 CONTINUE 
1000 FORMAT (1X.10F10.2) 

CALL fONTOR (MSZY»NnTF»FHlNc*LAB»* r NPT$L * FNPTSU*FRM*ZMtN,zMAX ) 

IXL = FLOAT 1 1 XL )+BLK+.5 

IXR=FL0AT(IXR)-BLK+.5 

CALL SETMIV ( I XL * I XR » I YB ♦ I YT ) 

IMMsO 

IMI=(MSZY-1 )*NOIF 
NFLAG1=1 
16 CONTINUF 

REWIND IRW 
READ (IRW) 

NPTSU=NPTSUU 

ntemp=nptsl 

nptsl=nptsu-npcre 

nptsu=nptsu+ ( NPTSU-NTEMP J-NPCRE 

18 CONTINUE 
STOP 
END 

SIRFTC GFTDA5 

SUrROUTiNp GFT5(DATA»NsPS»NSKiP,N(-Pw*NSCANO*iRw,rERR»NFLAG2* 
1NSTART»NRTLG*M0DF*NCRF*ITYPF.MSFC) 

DIMENSION DATAd ) *NDAT( 890) 

DATA NFLAG/O/ 

IF (NFLAG • NF • 0 ) GO TO 12 

NTEMP=(NCRE-1 )*NCPW 

NSCANOrO 

NL=36-NBTLG 

NRLNG=36/NBTLG 

nflagm 
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NSTMP=(NSPS+NSTART-1 )*NCPW+MSFC*19 

NSCANS*NSTMP 

NTMPsNSTMP/NBLNG 

n$tmp=nstmp/nblng*nrlng 

IF ( NSTMP *LT . NSC A NS ) NTMP=NTMP+1 
NSCANS=NTMP 

IF (NSKIP .EQ. 0) GO TO 4 O 2 
DO 401 I=1»NSKIP 
NSCANO=NSCANO+l 
CALL SKRBIN (IRW»1,RD ) 

401 CONTINUE 

402 CONTINUE 

12 CONTINUE 

IF CNFLAG2 .FQ. 0) GO TO 10 

13 CONTINUE 
M=(NSTART-1)*NcPW+MSFC*19 
IF (MODE .EG. 1 ) GO TO 50 

CALL REDTPR ( IRW*M0DE»IERR.NW»NSCANS.NDAT) 
GO TO 51 

50 CALL RFDTPc ( I RW »M0dE * I ERR *NW . NSCANS » NRAT ) 

51 CONTINUE 

1000 FORMAT ( 1 X ♦ 16 ) 

NNW=NW*6 
NFLAG2-0 
NSCAN0=NSCAN0+1 
10 CONTINUE 

DO 14 NN=I,NCPW 

IB=M+NBLNG 

TDsIR/NBLNG 

if=ir-nblng*td 

ibit=nbtlg*if 

MaM+1 

IF ( ITYPF .EQ. 1) GO TO 15 
ndat A=0 

CALL FLD(NDATA*NL.NpTLG*TRIT»NdaT ( ID) ) 
DATA(NN)=NDATA 
60 TO 14 
15 CONTINUE 
DATANsO.O 

CALL FLD(DATAN»NL»NRTLG»IBIT*NDAT ( ID) ) 
DATA<NN)=DATAN 

14 CONTINUE 
MaM+NTEMP 
RETURN 
END 


JIBFTC CTOR 

SUBROUTINE C0Nt0R(L»M,FHINC»LAB,FNPtSL,FNPT5U»FRM,2MIN,zMAX) 
DIMENSION I(5)*Vl(R)*V2(3)*Pl(5)»P2(3)»X(6)*YC5) 

DIMENSION V3(3) 

COMMON /Lab /Ni3»Nf 3»NT3»Z(2560) »KON(3400) ♦Ni2.NF2»NT2,un4O0) ♦ 
1 V(R400) 

LOGICAL JCODF 

IF (M.GT.O.AND.L.GT.O.AND.FHINC.GT.O.O.AND.LAB.GF.O) go TO 1 
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WRTTF 16*101) M,L*FHINC 

rfturn 

1 N=B400 

Nr= (ZMAX“2MIN ) /ARS ( FH T Nr ) 

FH=ZMIN 

2 CONTINUF 
NLAB=LAB-1 

IF (LAB .FQ. 0) NLAR=LAB+1 

DO 15 IH=1»NC 

Nl? = 6 

NF? = 5 

NT? = 5 

V<?)=FNPTSU 
U ( ? ) =fRm 
V(4)=FNPTSL 
U(4>=1.0 

K = 1 

NLAB=NLAB+1 

IF (NLAR .NF. LAB) GO TO 51 
K=? 

NLABrO 

51 DO 10 IL = 2 » L 
DO 10-IM=2,M 
1 (B) = < TL-1 ) *M+ I M 
1(?)=I(B)-1 

mwmi‘1 

i(5) = nn 

xn ) = im-i 

X { ? ) aIM-1 
xm = iM 
X(4)=IM 
X(5)*IM-1 
Y(1 )=IL-1 

Y ( ? ) = I L 

Y ( ■» ) = I L 

Y ( 4 ) = I L~1 
Y ( 5 ) a I L—l 
IZ?*I (2) 

121=1(1) 

!ZB=I (B) 

!Z4=T (4) 

TZ5=T (5 ) 

1000 FORMAT (1X*2F10,2*TM 
VI ( 1 ) = (X ( 1 )+X (4) ) /? , 0 
VI (2)=(Y(1 )+Y(?) )/2.0 

Vl (B)=(Z( IZ1)+Z( !Z2)+Z( IZ5)+Z( IZ4) )/4.0 
DO 4 J=l»4 
IZ=I ( J) 

V?(1 ) =X ( J ) 

V?(?)=Y{ J) 

V? ( B ) =Z ( T Z ) 

12=1 (J+l ) 

VB( 1 )=X( J+l ) 

VB(2)=Y( J+I ) 
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V8(8)=z< iz) 

IF (V 2(3) .If. 0.0 .AND. V8H) *LF. 0.0) GO TO A 
NV?=V?(8 ) + .5 
NV8=V8(8 )+.5 

TF t NV? .FQ. NV8) GO TO A 
1002 FORMAT ( IX *3F10, ? ) 

CALL STRiKF(V1»V2,V?*FH*P1 .P2*IFR) 

IF ( IFR .NF. 2 ) GO TO A 
IF (NF2+2 .LF. N) GO TO 59 
WRTTF (6*100) FH » N 
GO TO 60 

59 NF?=NF?+? 

U(NF?-1 )=PK?) 

VCNF2-1 )*Pl ( 1 ) 

ICON ( NF 2—1 ) =0 
U(NF?)=P? f 2 ) 

V(NF?)=P?(1) 

ICON { NF2 ) =1 

A CONTINUE 
10 CONTINUE 

60 IF (NI? .GF. NF?) GO TO 15 
JCODFs.TRUF. 

DO 61 J = N I ? *NF? 

jtfmp=j 

DO 61 JJ = ? * A ♦ 2 

IF <U(J) .FQ. U(JJ)) GO TO 68 
IP (V(J) .FQ. V(JJ)) GO TO 68 

61 CONTINUE 
NT?=NT2+? 

UtNT?-1 ) =U ( N I ? ) 

V(NT?-l)=V(Nl?) 

K0N(NT?~1 MO 

U ( NT ? ) =U ( N! ? + l ) 

V(NT?)=V(N1?+1 ) 

KON ( N I ? ) =1 
NI ?*N I ?+? 

GO TO 70 
68 J=JTFMP 

IF (KON(J) .FQ. 0) GO TO 62 
JJ=J~1 
JMAX=J 
GO TO 68 
6? JJ=J+1 
JMAX= J J 

68 VK1)=U(J) 

VI (?)*V( J) 

V? ( 1 )*U< Jj) 

V?(2)=V( JJ) 

JDIFsJMAX-Nl?-1 
DO 6A JSIJB=1*JDIF 
JDFX=JMAX-JSUB+1 

kdfx=jmax-jgub-i 

U(J0FX)=U(KDFX) 

64 V( JDFX)=V(KDFX) 

NT ?*N 1 2+? 
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IF (JCODF) GO TO 7? 

NT?=NT?+1 
U(NT?)=V? M ) 

VfNT?)=V?(?) 

KON ( NT2 ) =1 
GO TO 70 
72 NT?=NT?+? 

U(NT?-1 )=V1I1 ) 

V ( NT ?~1 )=V1 (?) 

KON(NT?-l ) =0 
U(NT?)=V?<1 ) 

V( NT? ) =V? ( ? ) 

KON ( NT? ) =1 

70 IF ( N 1 2 .GF. NF2) GO TO 16 
DO 71 J=N I ? * NF? 

JTFMP=J 

IF <U(J).NF.U(NT?),OR.V(J).NF.V(NT?) ) GO TO 71 

Jcodf=. false. 

GO TO 68 

71 CONTINUF 
GO TO 60 

16 CONTINUF 
NCNTssO 
J!=6 

DO 30 J=6*NT? 

JJ=J 

IF (KON(J) .FO. 1) GO TO 31 
IF (NCNT .GT. 4) GO TO 33 
DO 3? JI J=JI *JJ 
KON ( JI J ) =0 

32 CONTINUF 

33 NCNT=0 
J1=J 

GO TO 30 
31 NCNT=NCNT+1 
30 CONTINUF 

IF (K .NF. ?) GO TO 801 
CALL LABFL(FH) 

801 CONTINUF 

DO 800 J=6»NT2 
JUL=NXV(U( J) ) 

JVL=NYV(V(J)) 

JUR=NXV(U( J+l ) ) 

JVR=NYV(V( J+1 ) ) 

C WRiTF (6*1003) JUL » JVL * JUR ♦ JVR 

1003 FORMAT ( IX *416 ) 

IF ( KON ( J+l ) .NF. 0) CALL L I NFV ( JUL » JVL ♦ JUR * JVR ) 

800 CONTINUE 
15 FH=FH+FHI NC 

100 FORMAT (45H WARNING IN SURROUTINF cONTOR AT A HEIGHT OF * 
1F7.1 . 

246HNUMBFR OF POINTS FXCFFDFD ALLOWARLF NUM8FR OF ♦ 

316 */* 

421HFXFCUTI0N CONTINUING »//) 

101 FORMAT ( 28H FRROR IN SUBROUTING CONTOR * / 1 1 1 * 
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121H=NUMbfR Op grid rows ♦/.111* 

224H=NUMrfR Op GRID COLUMNS ,/,F11.3» 
333H=DISTaNCF RFTWFpN CONTOUR LARFLS ♦/♦111. 
429H=FRFQUFNCY OF CONTOUR LARFLS ♦/♦ 
521HFXFCUTTON TfRMINaTFO .// ) 

RETURN 

END 


SIRFTC STRK 

SUPROUT tNp STR JKF(V1 .V2.V3.FH.Pl.P2. TFR > 

DIMENSION Vl ( 1 ) »V2( 1 ) »V3(1) ,Pl ( 1 ) »P2( 1 ) »Q(3.3) 

DO 6 1=1 .3 
Pi ( I)=0.0 
P?U ) = 0, 0 
Q(! .1 )=V1 ( I ) 

0(1 * 2 ) =V? ( I ) 

6 0(1 ,3)*V3 ( I ) 

1000 FORMAT (1X.3F10.5) 

DO 7 1=1 *3 
II = ! 

DO 1 J=! I * 3 

IF (Q(3*TI)“0(3.J)) 1.1,2 
2 DO R L=1 ,3 

A=0(L.II) 

0(L,II)=Q(L.J) 

0(L, J)=A 
R CONTINUF 

1 CONTINUF 

7 CONTINUE 
I FR = 0 

IFCOC3.3j.LT.FH.OR. Q ( 3 , 1 > . 6T . FH ) RETURN 

IF (0(3.1 ) .FQ.FH.AND.O(3,2) .FQ.FH. AND.Q( 3.3) .EQ.FH) RETURN 
As ( PH~0( 3*1 11/ (0(3*1 ) -0(3.3) ) 

Pi (1 )=(Q(1 ,1 )-Q(l .3 ) >*a+Q( 1 *1 ) 

Pl(2) = (0(2.1) -0(2.3)) *A+0( 2 *1 ) 

P) ( 3 ) =FH 

IF (0(2.1 ) .NF. FH) GO TO 3 
IF (0(3,2) .NF. FH) GO TO 4 

P? ( 1 ) =0 ( 1 » 2 ) 

P? ( 2 )=0( ? . 7 ) 

P?(3)sO( P.2) 

IFR = ? 

RETURN 

3 IF (0(3,2) ..GT. PH) GO TO 5 
A= ( FH-0 (3.2) )/(Q(3.2)-Q(3.3) ) 

P2(l ) = (Q( 1 ,2)-0(l ,3) ) *A+Q (1.2) 

P? ( 2 ) * l 0 ( 2 ♦ 2 ) -Q ( 2 ♦ 3 ) )*A+0(2.2) 

P?(3)=FH 
IFR*2 
RETURN 
4 IFR=1 

P? C t J =P 1 (1 ) 

P 2 ( 2 ) =P 1 (2) 

P?(3)=FH 
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RFTURN 

5 A=(FH”Q(3*1) )/<Q(3,i )-Q(3*2) ) 

p?m = iQ(i»ii-on»?i ) * a+q ( i , i ) 

P?{2)s(OC?>])-Q(7»?) )*a+Q(2»1 ) 

p?n)=FH 

IFR = 2 

RETURN 

END 


JIBFTC MARK 

SURROUTINF LABFLIH) 

COMMON /LAR /Nt 3 »Nf 3 *NT 3 ,Z (2 560 ) »KON ( 3 A 00 ) * N I 2 . Nf 2 ^ 2 *U( 3400 ) ♦ 
1 V(BAOO) 

1 2 = 5 

12 11=1241 

DIS= 0.0 
15*1141 

IF (1144 .GT. NT 2 ) GO TO 99 

DO 10 T 2 =IS*NT 2 

IF (K 0 NU 2 ) .FQ. 0 ) GO TO U 

10 DlS=OTS+SORT( (U( T 2 )-U( 12 - 1 )) ** 24 t V ( t 2 ) -V ( ! 2 - 1 ) )** 2 ) 

11 ! 1 2 = 1 2-1 

IF (DIS .LT. 270 . 0 ) GO TO 12 
lM=UI 24 ll )/2 
I r» T F= ( 1 1 2- 1 1 ) /? 

IF ( ID! F .LT. ?) GO TO 12 
DO ?0 J= 1 » IDIF 
I L= I M — 1 
I U* I M-f J 

OTS*SQRT( (UnU)-U(TL) ) **24 ( V ( I U ) -V ( I L ) )**2) 

IF (DIS .GT. 68 . 0 ) GO TO 16 
20 CONTINUF 
GO TO 12 

16 I 6 =IL 41 

DO 19 JJ=IS»IU 
19 KON ( JJ ) sO 

RFALM= 10 . 0 F 410 

IF (U(IU)-U(lL) .Nf. 0.0) RfaLM=(V( TU)“V( IL) )/(U(IU)-U( lU) ) 
ANGLsATaN(RFALM) 

R=(DTS- 68 . 0 )/ 2.0 

Ss6P.04R 

IF CUt I U 1 .GT. U ( IL ) ) GO TO 30 
1 ! = IU 
IU*IL 
I L= 1 I 

30 XP=(R*U( TU)+S*U( IL) ) /DIS 46 . 8 *Sl N ( aNgL ) 

YP=(R*V( IU)+S*V( IL) ) /DIS- 6 . 8 # C 0 S( ANgL) 

DO 18 JJ= 1*5 
ND= 5 -JJ 

IF (ARS(H) .LT. 1 0 . 0 **( JJ- 1 ) ) GO TO 17 
18 CONTINUF 

17 CONTINUF 
IXPsXP 
IYPsYP 
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CALL LABLVCH* IXP*IYP*ND»1 .5) 
GO TO 1? 

99 RFTURN 
END 


SORTGIN ALPHA»SYSUT2»RPW 

SIRFTC MOD6 

subroutine bwndRs 

DIMENSION X( 12*256) *Y ( 12*256) »MfHAN( 12) »NN( 256) *KSYM(49) » JSYM( 
1256) 

DIMENSION NWH I EH (12) 

NAMELTST/lNPUT6/NSrANS*NSTART»NSPS*NrH* N VAR*NSYM*iSUM,NRTLG* 
1M0DE*ITYPF.MSFC .NSKiP,NbLK,iNcX.INcY.NSTX,NSTY,NcRf 
NaMELIST/NCHUSE/NWhTCH 
EQUIVALENCE (NSCAN,NSCANS) 

EQUIVALENCE (NSTRT,NSTART) 

EQUIVALENCE (NCOL,NSPS) 

EQUIVALENCE (NCHAN.NCH) 

ICARn=5 
IPR 1 NT=6 
INTAPE=10 
I0TAPE = U 

READUCARD»INPUT6) 

WRITE ( I PR I NT * I NPUT 6 ) 

REAO(ICARD*NcHUSF) 

WRITE(IPRINT*NCHUSE) 

1 FORMAT ( 1 X * 7 1 4 ) 

8 F0RMAT(1X*12H ) 

READ! I CARD *5 ) (KSYM( t ) *1 = 1 » NSYM ) 

5 FORMAT ( 1 X » 60a 1 ) 

NELAG=0 
AVF= I SUM 
APOPaO.O 
DXAVFsO • 0 
DYAVF*0.n 

dzave=o.o 

NSAVsNSCAN 

IF CNSKlP *EQ. 0) GO TO 98 

DO 97 1 = 1 * NSK. I P 

CALL SKRBIN( INTAPE*! *NOP) 

97 CONTINUE 

98 CONTINUE 

2 FORMAT ( 1 HI ) 

A FORMAT CSX *111 II ) 

11 = 1 

KK=NSTRT-1 

160 IFC II.EQ.NSCAN) GO TO 510 
11=11+1 
NFLAG2=1 
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IFUI.NF.2) GO TO 290 
00 170 JJ=1 *NCOL 

CALL GET 6 ( X ( 1 * JU) »NrOL#0*NcHAN»NSCANO,iNTAPF*IFRR*NFLAG2*NSTRT» 
lNRTLG«MOnE*NcRF* TTYPF»MSFC) 

170 CONTINUE 
290 CONTINUE 
NFL AG?=1 

DO 300 JJ=l,NCOL 

CALL GET 6 ( Y ( 1 »JJ) »NrOL*0*NrHAN,NSCANO* iNjAPr* IFRRtNFLAG2»NSTRT» 
1NRTLG*M0DF*NCRP» iTYPEtMSFC) 

300 CONTINUE 

DO 380 JJ=2»NC0L 
I J= JJ-1 
XSUM=0.0 
YGUMsO.O 
ZSUMaO.O 

DO 360 TCHAN=1,ISUM 
1 1 CHN=NWHI CH ( I CHAN ) 

XDIFF=Y( I ICHN*JJ)-Y(1ICHN,I J) 

YDIFF=Y( IICHN» JJ)-X( ITCHN,JJ) 

xsum=xsum+xdiff*xdiff 

ysum=ysum+ydiff*ydtff 

zsum=zsum+xdiff*ydtff 

360 CONTINUE 

XSUMsXSUM/AVF 

ysum=ysum/ave 

ZSUMaZSUM/AVF 
APOP=APOP+1.0 
AA=1.0/APOP 
BB*1 . O-AA 

dxavf=rr*dxavf+aa*y.sum 
OYAVc=RR*DYAVF+A A*XSUM 

OZAVF=BR*DZAVF+AA*ZSUM 
XSUM*S3RT (XSUM^ 

TF (XSUM .LT. 53.4) GO TO 365 
364 XSUM=53 . 0 
YSUM=53 • 0 
GO TO 366 
364 CONTINUE 

YSUMsSORT (YSUM) 

IF (YSUM .LT. 53.4) GO TO 366 
GO TO 364 
366 CONTINUE 

CALL JNTPB(YSUM,XSUM,NFLAG»0*0*KSYM{3) ,NPOP» 

2DXAVE.DYAVE* 

3D7AVF. 

I JJ » NSCAN » 1 1 »NC0L ♦ JSYM.X »NSTRT»NN, INCX ) 

N5CAN=NSAV 
380 CONTINUF 

DO 500 JJ=1,NC0L 
DO 490 ICHAN = 1 * I SUM 
1 1 CHN = NWHI CH ( ICHAN ) 

X C I I CHN ♦ JJ ) =Y ( II CHN ♦ J J) 

490 CONTINUE 
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500 CONTINUE 
GO TO 150 
510 CONTINUE 
NFL AG* 1 

CALL JNTPB<YSUM,XSUM,NFLAG.0»0»KSYM(3) ,NP0P, 
2DXAVF *DYAVE * 

3DZAVF, 

1 JJ*NScaN » 1 1 *NCOL» J$YM,X»NSTRT»NN» INCX) 

RFWIND IOTAPF 
REWIND TNTAPF 
C CALL CLEAN 

RETURN 
END 


$IRFTC PLOTT6 

SUrROUTtNF LARFL6(NSTART»NSTOP* incRf 1 
DIMENSION IOUT(1?0) 

NdTF=(NSTOP-NSTART+1 1/INCRF 
11=0 

DO 1 I=NSTART.NSTOPtlNcRE 
II=I 1+1 

IOUTUDsT/lOOO 
1 CONTINUE 

WRtTf ( 6 * 1 0 } (T0UT(T)»I=1.NDTF1 
11=0 

DO 2 I=NSTART»NSTOP,lNCRE 
11 = 11+1 

1 OUT (I I 1=1/100-1/1000*10 
? CONTINUE 

WRtTf ( 6 ♦ 1 0 ) ( I OUT (I) *1=1* Nd IF) 

11=0 

DO 3 I=NSTART ,NSTOP,lNcRE 
IT=I 1+1 

IOUTUI ) = I/10-T/100*10 

3 CONTINUF 

WRITE ( 6 « 1 0 ) ( I OUT (I )*! = !» ND IF) 

11=0 

DO 4 !=NSTART»NST0P,INCRE 
11=11+1 

IOUTUI 1 = 1-1/10*10 

IF (IOUTUI) .LF. 0 ) TOUT ( T I ) = 0 

4 CONTINUF 

WRTTF ( 6 » 1 0 ) (TOUT U 1 , 1 = 1 *NDIF 1 
10 FORMAT (11X»l?0ll) 

RETURN 

END 


SIRFTC GFTnA6 

SUBROUTINE GFT6( DAT A » NsP$ * NSK I P * NfP » NSC aNO * lRw» I ERR *NfLaG 2 ♦ 
lNSTART»NPTLG»MODF»NrRr , lTYPr,MSFC 1 
DIMENSION DATAU)*NDAT( 890) 

DATA NFLAG/O/ 

IF (NFLAG ,Nf. 0 ) GO TO 1? 
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NSCANOaO 

NTPMP= (NrRF- 1 )*NCPW 
NL = 3f>-NBTLG 
NBLNG=36/NBTLG 
NFLAG=1 

NSTMP=(NSPS+NST ART-i )*NCPW+MSFC*19 

NSCANS=NSTMP 

NTMP=N$TMP/NBLNG 

NSTMP=NSTMP/NBING*NRLNG 

IF ( NSTMP , LT . NSC A NS ) NTMPsNTMP + i 

nscans=ntmp 

IF (NSKIP .FO. 0) GO TO 4 O 2 
DO 401 I=l»NSKlP 
NSCANO=NSCANO+l 
CALL SKRBIN { IRW, 1 »RD > 

401 CONTI NUE 
40? CONTI NUF 

1 ? continue 

IF (NFLAG2 .FQ. 0 ) C,0 TO 10 

13 CONTINUE 

M={NSTART- 1 )*NcPW+MSFC*19 
IF (MODE .EG. 1 ) GO TO 50 

CALL RFDTPR ( IRW.MODE, IFRR,NW,NSCAN5*NDAT) 
GO TO 51 

50 CALL RFHTPC (TRW»MOnF,!FRR*NW»NSCANS»NOAT) 

51 CONTINUE 

1000 FORMAT ( 1 X * 1 6 ) 

NNW=NW*6 
NFL AG?=0 
NSCAN0=NSCAN0+1 
10 CONTINUE 

DO 14 NN= 1 *NCPW 
!BsM+NBLNG 
id=ib/nblng 
if=tb-nblng*id 

IBIT=NBTLG*IE 

M=M +1 

IF UTYPF .FQ. 1) GO TO 15 
NDAT a=o 

fall FLn(NOATA*NL*NRTLG»TBTT*Nr)AT(TD)) 
DATA(NN)=NDATA 
GO TO 14 
15 CONTINUE 
DATAN=0. 0 

CALL FLD(DATAN,NL»NbTLG. IBIT*NDAT( ID 1 ) 
DATA(NN)=DATAN 

14 CONTINUE 
MaM+NTEMP 
RFTURN 
FND 


SIBFTC JOINT 

SUBROUTINE JNTPB(DATAH*DATAV»NFLAG*MIX.Miy»ALPNuM,NPOP, 
2DXAVF.DYAVF* 
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3DZAVF« 

1JJ,NScaN* TSCAN,NrOL»lSYM,Nx*NSTRT*NN, INCXY) 

DIMENSION INCXY(l) 

DIMENSION NP(54»54) 

DIMENSION DATAI1?) 

DIMENSION IRI N ( ?55 I 
DIMENSION ISYM(i ) ,tOUlT( TOO) 

DIMENSION Nxn ) ♦NNO ) 

DIMENSION ALPNUM< 1 ) ♦ALPHA(1?0I ,rORDX(3) 

DOUBLE PRECISION A ( ? * ? > *E I OFN ( 2 » 2 ) 

INTEGER ALPNUM,ALPHA»BLANK 

data ASTRIK/ih*/ 

DATA XMARK/lHX/ 

DATA BLANK/6H / 

DATA NFLAG4/0/ 

1060 FORMAT {48X»25HDATA SWITCH HAS OrcURRFD ) 

1061 FORMaT(APX»30hJOINT PROBABILITY DISTRIBUTION ) 

106? F0RMATMH1 ) 

1063 FORMAT (44X»11h X-AXtS IS »i6«6X«11h Y-AXlS IS »2l6) 

1066 FORMAT ( 3 OX *6HDX AVF= * Fl 5 . 7 , 6HDYAVF= , Fl 5 . 7 , 6H0ZaVF= , F l 5.7 ) 

1040 FORMAT ( 1 H » 67H MAXIMUM PROBABILITY OF UNcOMMONaL J T Y EXCEEDED" rONT 
1INUE EXCUTION *216) 

1064 FORMAT UX*26HSYMBOL N/SYMBOL ) 

1065 FORMaT(11X,121(1h*)*/*11X.1h*,55X*5hPaRT *il,4H Of * 1 1 , 53X , 1H* . / 

1 *lix*121 < lH*n 

IF I NFL AG 4 ,GT. 0) GO TO 80 
NFLG3=0 

1000 F0RMAT(1X»47A1 ) 

NFLG=0 

NI=? 

I RW= 1 
NFLGFN=0 
IRT=11 
DO 1 1=1 * 54 

DO 1 J= 1 * 64 
NP ( I * J ) =0 
1 CONTINUE 

REWIND IRT 
REWIND IRW 
NFLAG4=1 
80 CONTINUE 

IF INFLAG. GT. 0) GO TO 13 
NC=0ATAV+1 . 5 
NR=DATAH+1.5 
IF (NC .LT. 1 ) NC=1 
IF <NR .LT. 1 ) NR= i 
NP { NR *NC ) =NP { NR » NC ) +1 
1=84*1 NC“1 ) +NR 
1BTN( JJ)=I 

IF (JJ .LT. NCOL) GO TO 18 
tBTNm = IRlN(? ) 

WRITFI IRW) I IRINI I I) *1 1 = 1 *NCOL) 

15 CONTINUE 
RETURN 

13 CONTINUE 
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17 CONTINUE 

REWIND IPW 
I OPT *» 1 

I Ns ? 

I Me ? 

RhO=1 .0/(10, 0**5) 

AM ,1 ) =DXAVE 
A ( ? ♦ 2 ) =DY AVE 
AM *7)=DZAVF 
A(?#1 )=DZAVF 

CALL DJCOBI ( A*!M*iN*IOPT,RhO»eRR»ETGFN) 

C WR I TF ( 6 » 1 067 ) A(l»l) *A(1*2)*PTGFN(1*1)*EICFN(1,2) 

C WR I Tf ( 6 * 1 067 ) A(2»1)*A<2.2)#EIGFN(2#1).FIGFN(2»2) 

1067 FORMAT ( 1X*2E15.7*10X*2F15.7) 

DXAVF = 1.0/AM #1 ) 

DYAVF=1.0/A(?»?> 

A ( 1 ♦ 1 )=F1GENM »1 )*DXAVf 
AM *2)=FIGFN(2*I )*r»XAVF 

A(2*l )=E1GEN(1,2)*DYAVF 

A(2»2)=FlGFN(2t?)*DYAVF 

DXaVE=EIGEN(1*1 )*A( 1 ♦1)+EIGFN(1.2)*A(2*1 ) 

DYaVE = ETGFN (2 * 1 ) *A ( 1 *2)+FlGFN (2»2)*A(2#2) 

DZaVf=ETGFN(1#1)*A(1*2)+EI gfN 1 1 *2 ) * A ( 2 »2 ) +FIGFN < 2 * 1 ) *A ( 1 » 1 > 
1+FTGfN(?»?)*A ( 2 » 1 ) 

WR T TF ( 6 ♦ 1 066 ) OX AVf,OYAVF *DZAVF 
1=0 

DO 150 NC = 1*64 
DO 150 NR = 1 * 54 
1 = 1 + 1 

IF (NP(NR.NC) .EQ. 0) GO TO 130 
XXX =NR*NR 

yyy=nc*nc 
ZZZ =NR*NC 

SUM = nXAVF*XXX+DY A Vf*YYY 
1+DZ AVE*ZZZ 

I F ( SUM.GF ,1.0) GO TO 115 
NX ( I ) =0 
GO TO 130 
116 NX ( I ) =~1 
130 CONTINUE 

WRlTF ( 6 » 1 06? ) 

WR T TF < 6 » 1 061 ) 

WR T Tc ( 6 * 1 066 ) DX aVf *DYaVF #DZ AVf 
WR lTFI6»1 064) 

C CALCULATE TABLE 

MAXKNT=0 
DO ?? NC=1 #53 
DO 2? NR=1 #63 

IF ( NP ( NR #Nr ) ,GT. MAXKNT) MaXKNT=NP ( NR ♦ NC } 

2? CONTINUF 

tF (MAXKNT .LT. 46 ) MAXKNT=46 
NFACT=MaXKNT/46 
XXX=FLOAT (MaXKNT ) /46.0 
NFAC=0 

IF (NFACT .LT. 1) NFACT=1 
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WRTTF(6.l050) bLANK.NFAC 
NFAC=NFAC +NFACT 
f>0 23 1 = 1 *46 

WR|TF(6»l050) aLPNUM! I ) »NFAC 
1050 FORMAT (3X.A6.6X. 16) 

nfac=nfac+nfact 

23 rONTTNUF 

WRITF(6*106?) 

C PRINT DISTRIBUTION ON PAGE 

CORDXU 1=0.0 
C0RDX(2)=r0R0X( 1 )+60.0 
COR OX ( 3 ) =CORDX (1 ) +11 0.0 
1021 FORMAT C 1 HI ) 

DO 65 I £ND= 1 * 54 
NC=56-IEND 
DO 66 1=1 .54 
ALPHA ( I ) =RLaNK 
66 CONTINUE 

IF (NFLG .EQ. 1 ) GO TO 69 

DO 68 1=1.54 

TB1N1I)=0 

68 CONTINUE 

69 CONTINUE 

DO 64 NR = 1 ,54 

XX a FLO AT ( NP ( NR »NC ) ) 

IF (NFLG .EO. 1 ) GO TO 91 
IBIN(NR ) =NP (NR »NC ) 

91 CONTINUE 

lCAR=XX+( 1 .001-1 .O/XXX) 

IF (ICAR ,GT. 46 ) ICAR=46 
ALPHA (NR )=ALPNUM{ ICAR) 

64 CONTINUE 
COROY=FLOAT(NC ) 

YMARG=XMARK 

IF (NC .NE. 54-IENO/10*l0) YMARG=ASTR IK 
WRITF(6»1008) cORDY.YMaRG.I ALPHA! I 1.1=1*54) 
1008 FORMAT! IX » F8. 1.2X » A1 .120A1 1 

65 CONTINUE 
WRITE (6.1010) 

WRtTf( 6*1 01 1 ) rORnxn ) *cORr>X ( 2 ) »cORoX ( 3 ) 
1011 FORMAT (6X »F10.4.50X»F10,4,40X.F10.4) 

nflg=i 

WRlTF(6.1 062) 

N$UB=1 
LWFR=1 
LOW=NSTRT 
705 CONTINUE 

NHI =L0W+1 20-1 

nuppfr=lwfr+i 20-1 

IF (NUPPFR ,gt, NCOL) NUPP er=ncol 

WRTTF(6.1 062 ) 

CALL LABFL6(L0W,NHI *1 ) 

DO 131 1 1 =Nl .NSCAN 

RFAD(IRW) ( IRIN! JJ) , JJ=l .NCOL) 

DO 135 J J= 1 .NCOL 
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irHFrK=TRTN( JJ) 

Jchfck=nx(ichfcK) 

! F (JCHFCK .NF. 0) r,0 TO 117 
ISYM( JJ)=ALPNUM(NStJR-1 ) 

NN( JJJ.O 
GO TO 135 

117 I SYM ( JJ ) =ALPNUM ( NSUB ) 

NN(JJ)=-1 
136 CONTINUE 

IF (NFLGFN .NF. 0) GO TO 116 
WRITF(IRT) INN( JJ) ♦ JJ=1 ,NCOL ) 

CALL PLTrf 6< lSYM,NfOL,NeLK, jNcXYt 1 ) * lNcXY( 2) *lNrXY( 3 ) , 
1INCXYU) *NCRF) 

136 CONTINUE 

WRlTF(6*l0i6) II*(ISYM(JJ) , J J=LWFR * NUPPFR ) 

1036 FORMaT(5X»I6*1?OA1 ) 

1036 FORMaTMX* 16) 

131 CONTlNUF 

nflgfn=i 

REWIND IRW 

LWFR=NUPPER+1 

L0W=NHI+1 

IF (NUPPFR .LT. NCOL) GO TO 705 
905 CONTlNUF 

lOlO FORMAT (iiX»1?(10hX*********)) 

nflgfn=o 

return 

END 


SJBFTC PLLT6 

SUBROUTINE PLTBF6< I 5YM , NN , NsOR * INCX * I Nc Y * N$TX » NsTY , NCRF ) 
DIMENSION ISYM( 1 ) jNBUFFR ( 50 ) 

data nflag/o/ 

DATA NFLAGl/O/ 

IF ( NFL AG .NF, 0) GO TO 10 
CALL CAMRAV ( 35 ) 

CALL BUTTV(l) 

10 CONTlNUF 

IF (NFLAG1 .NE. 0) GO TO 20 
CALL FRAMEV(O) 

NC0UNT=0 

INCRX*NSTX 

incRy=nsty 

nflagi=i 

IF (NFLAG .EQ. 0} GO TO 11 

call aPRNTV ( INcX» INCY*NN,NBUFER» INcRX*INcRY) 

INCRY=INCRY+INCX 

INCRX=INcRX+IARS( TNGY ) 

NC0UNT=NC0UNT+1 

11 nflag*i 

20 CONTINUE 

ncount=ncount+i 

DO 1 1=1*50 
NBUFFRI I ) =0 
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1 


CONTINUE 
I A*0 

DO ? 1*1 *NN 
lBstA+6 
I C= I B/6 
!D=IB-IC*6 
IE=ID«6 
I A=1 A+l 

CALL FLdCNRUFFRC IC) *TF» 6. 0*ISYMCI>) 

2 CONTINUE 

CALL APRNTVnN<"X*lNrY*NN,NBUFFR»lNCRX»lNcRY> 
incRy=incry+incx 

1NCRX = INCRX+IABSUNCY) 

IF (INCRX .GE. 1023) NFLAG 1=0 

RETURN 

END 


SORIGIN ALPHA »SYSUT2*RFW 

SIBFTC MOD? 

surRoutine clasfy 

COMMON /LAR1/XraR(A2»12) .SIGMA C A2 » 12 > * ROT { 42 » 12 » 1 2 ) 

COMMON /LAR2/X( 12) »ALPhA( 49) *NsPS »NSCANS ♦ NCHAN , LT9 , LtIP *LT 11 *LT12* 
1LT13»LT1 * IXXX»IYYY, 

INSTART * NS TOP* 

inrtlg»modf*i type* msec. I4*NCRE* 

1 NSKIP»INcX» INcY*NSTXiNSTY 
NaMELIST/INPuT7/NPaSS.NcLUST 

Namelist/ INPUT a/NSPS»NScaNS,NcH* LTl.LT9»LT10.LTll*LT12»LT13 t 
lNSTARTtNSTOP,NBTLG«MOr)E* ITYPE*MSFC» I A»NCRE,N5KIP,INCX,INCY.NSTX* 
2 NSTY*IXXX» IYYY 
EQUIVALENCE (NCH.NCHAN) 

Read c s » i NPUT 7 ) 

RFAD ( 5 * I NPUT A ) 

WRlTF( 6 »TNPUTA) 

READ ( 5 ♦ 1 006 ) ( ALPHA t I ) * I =1 *48 ) 

1006 FORMAT ( IX * 60A1 ) 

kount=nclust 

nscans=nscans-i 

IN1TCL=NCLUST+1 

DO 1 1=1 *NPASS 

CALL TRUCK(NcLUST»NPASS ) 

CALL SFQMRG (NCLUST»KOUNT*INlTCL ) 

CALL CLASS (KOUNT , I * NPaSS ) 

NCLUST =KOUNT 
IN1TCL=K0UNT+1 
1 CONTINUE 
RETURN 
END 
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*IRFTC GFTDA7 

SUpROUTlNE GET7(DATA»NsPS»NSKlP,NcPw»NSCANO, iRw, IERR»N fLaG 2. 
lNSTART»NBTLG»MOr)E,NrRF*lTYPF,MSFC ) 

DIMENSION DATA ( 1 ) »NDAT ( 890) 

DATA NFIAG/O/ 

IF (NFL AG .NP. 0 ) GO TO 1? 

NSCANOaO 

NTFMP=NCPW*(NCRE-1 > 

NLa36~N8TLG 

NBLNG=36/NBTLG 

NFLAG=1 

NSTMP=(NSPS+NSTART-i ) *NCPw+MsFC* 1 9 
NSCANS=NSTMP 

ntmp=nstmp/nblng 

NSTMP=NSTMP/NBLNG*NRLNG 

IF ( NSTMP .LT. NSCANS) NTMP=NTMP+1 

nscans=ntmp 

IF (NSKIP ,EQ. 0) GO TO 40? 

DO 401 1=1,NSKIP 

NSCAN0=NSCAN0+1 

CALL SKRriN ( I RW * 1 ,RD ) 

401 CONTINUF 

402 CONTINUE 
1? CONTINUF 

IF (NFL AG 2 .EQ. 0) GO TO 10 
13 CONTINUE 

M= (NST ART-1 )*NcPW+MSFC*19 
IF (MODE .EQ. 1) GO TO 50 

CALL REDTPR ( !RW»MODE*IFRR*NW.NSCANS*NDAT ) 

GO TO 51 

50 CALL REOTPc ( iRW *M0 dE ♦ I ERR ♦ NW.NSCANS * NpAT ) 

51 CONTINUE 

1000 FORMAT (IX. 16) 

NNW=NW«6 

NFLAG?sO 

NSCANO=NSCANO+) 

10 CONTINUE 

DO 14 NN=l . NcPW 

IBsM+NBLNG 

ID=IR/NBLNG 

IF=IB-NBLNG*ID 

IBITsNBTLG*IE 

M*M+1 

IF ( ITYPF ,FQ. 1 ) GO TO 15 
NDATA=0 

CALL FLD(NDATA»NL»NRTLG»IBIT»NDAT( ID) ) 

DATA(NN)=NDATA 
GO TO 14 
15 CONTINUE 
DATANsO.O 

CALL FLd( DATaN.NL »NrTLG* IB IT *NDAT( ID) ) 

DATA(NN)=DATAN 
14 CONTINUE 
MsM+NTEMP 
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return 

END 


SIBFTF PL0TT7 

SURROUTiNf LARC17 (NSTaRT »NSTOP*INcRf ) 

DIMENSION IOUT(1?0] 

NDTF=(NSTOP-NSTART-t.l )/INcRf 
ii=o 

DO 1 I=NSTART»NSTOP .INCRE 
11=11+1 

TOUT C 1 1 ) =1 /I 000 

1 CONTINUE 

WRlTF (6*10) ( f OUT (T)*I=1*NDTF) 

11=0 

DO 2 1=NSTaRT,NST0P*INcRE 
11 = 11+1 

10UTUT) =1/100-1/1000*10 

2 CONTINUE 

WRITE (6 * 1 0 ) (I0UT(I)*! = 1*N0TF) 

11=0 

DO R I=NSTART*NSTOP,INcRE 
11 = 11+1 

!OUT(II)*I/10-I/100*10 
R CONTINUE 

WRITE (6*10) ( IOUT ( I ) * 1 = 1 *NDIF) 

11=0 

DO A I=NSTART*NSTOP,lNCRE 
1 1 = I I + 1 

!OUT(II)=I“I/lO*lO 
IF ( IOUT (II) .LE. 0 ) I OUT (II ) =0 
A CONTINUE 

WRITE (6*10) (I0UT(I)*I=1»NdtF) 

10 FORMAT MIX* 12 Oil) 

RETURN 

END 


SIRFTC PLLTT7 

surRoutine pltrf7< isym,nn,nsor* incx*incy,nstx.nsty*ncRe* 
1 nflag»nflagi ) 

DIMENSION ISYM( 1 ) ,NrUFER( 50) 

DATA NFLG/0/ 

IF (NFLG *NE* 0) GO TO 10 
NFLG=1 

CALL CAMRAV(R5) 

CALL RUTTV ( 1 ) 

NSOR=102A/T ARS( INCY+INCX ) 

10 CONTINUE 

IF (NFLAG1 .NE . 0) GO TO 20 
CALL FRAMEV(O) 

NCOUNT=0 

incrx*nstx 
l ncRy=nsty 
NFLAG1*! 
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11 

20 


1 


2 


IF ( NFL AG .EQ* 0) GO TO 11 

CALL aPRNTV ( iNcX*!NCY*NN,NbUFPR 

INCRY=INCRY+INCX 

INCRX^INCRX+IABSHNCY) 

NCOUNT=NCOUNT+l 


♦INcRX*INcRY) 


NFLAG=1 
CONTINUE 
NCOUNT cNCOUNT + 1 
DO 1 1=1*50 
NBUFERd ) =° 
CONTINUE 


I A = 0 

DO 2 1 = 1. *NN *NCRE 
IB = I A+6 
IC=IB/6 
T D* I R-IC*6 
IE=ID*6 


IA=1 A+l 

CALL FLD(NBUFFRUC1 *IE*6* 0*ISYM(I») 

CONTINUE M „ 

CALL APRNTV ( lNcX*!NrY»NN,NRUFER»lNCRX»lNcRY) 
I NcRY=I NCRY+I NCX 
1NCRX=INCRX+TARS( TNrY) 

IF (INCRX .GE. 1024) NFLAG1=0 

RETURN 

END 


JORIGIN BRAVO,SYSUT2*REW 

$ I RFT C M0D7A 

SURROUTINF TRUfKtNrrNT.NPASS ) 

DIMENSION NN ACC (12«?56)»MTA8fll)*lPRT(256) * I PLOT ( 256) 

DIMENSION NT8L(400) 

COMMON /LABl/XRARU2*12) .SIGMA (42*12 )*ROT 142*12* 12 ) 

COMMON /LAB2/X112) »NSYM (49)»NSPS»NSCANS*NchaN»LT9*LT10»LT11*LT12 
1LT13»LT1 , IXXX,IYYY, 

1NSTART*NST0P, 

inrtlg»mode*itype.msfc*ta»ncRe» 

1NSKIP»INcX* INCY*NSTX*NSTY 
NFLGXX=0 
NFLAGX=0 
REWIND LT11 
REWIND LT1 
NFLAG1 =0 

mftn=o 


IXIY=1XXX*IYYY 
DO 10 1=1*1 YYY 
MT AB ( I ) = T 
10 CONTINUE 

DO 50 1=1 *400 
NTBLU ) = 1 
50 CONTINUE 

DO 11 1=1 * I YYY 

RFAD ( LT 1 1 ) (NNACCd *JJ) »J J = 1 *NSPS) 
mfin=mfin+i 


f 
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11 CONTINUF 

NUP=NSPS-IXXX+1 
NCNT=NCCNT+l 
NFLAG*0 
200 CONTINUE 

in=MTAB<l) 

DO no Jj = l »NSPS 

IF (JJ .GT • NUP) GO TO 102 

I J=JJ 

JlsJJ+IXXX-1 
NZERO=0 
IKNTsO 
I SUM=0 
JI JeMTAR ( 1 ) 

NTFMP= NCCNT+1 
DO 101 I=IJ»JI 
DO 100 JI J«1 , IYYY 
1!J=MTAB(JIJ) 

IF (NNACCUIJ.I) .LF. NCCNT .AND. NNACC ( I I J • I ) . NE. 0) GO TO 102 
IF (NNACCUlJtin 102*107.106 

106 IF (NNACCHU.n .GT. NTFMP) NTeMP=NNaCC ( 1 1 J* I ) 

GO TO 100 

107 NZFR0=NZFR0+1 

100 CONTINUE 

101 continue 

IF (NZERC .NE. I X I Y ) GO TO 1 O 5 
DO 103 I=IJ.JI 
DO 104 JI J = 1 * I YYY 
NNACCUU.I )=NcNT 

104 CONTINUE 
103 CONTINUF 

NCNT=NCNT+1 

IF (NCNT ,GT. 400) GO TO 999 
GO TO HO 

105 CONTINUE 

DO 108 l = U*JI 
DO 108 JI J=1 * I YYY 

IF (NNACC( JI J. I ) .FG. 0) NNACC ( JIJ. I) =NTEMP 

108 CONTINUF 
GO TO HO 

102 CONTINUF 

110 continue 

DO 111 JJ=1 >N5PS 

IF (JJ .EQ. NSPS ) GO TO m 

IF ( NNACC ( I 1 1 * JJ ) .LE. NCCNT) GO TO 111 

IF (NNACC(ITI.JJ+1) .LF. NCCNT) GO TO 111 

IF ( NNACC fill* J J 1 , LF. 0) GO TO 111 

IF (NNACC C I I I ♦ JJ+1 ) .LF. 0) GO TO 111 

IF (NNACC( I 1 1 . JJ) .EQ. NNACC ( I I I » JJ+1 ) ) GO TO 111 

IJsNNACCUI I.JJ) 

J!=NNACC(III»JJ+1) 

IF (JI .GT. 400 .OR, IJ ,GT, 400) GO TO 111 
IF (NTRL(Jl) .GT. NTRLUJn GO TO 125 
NTRL ( I J ) aNTBL ( Jl ) 

GO TO 111 
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125 CONTINUE 
NTRLCJI ) =NTrL ( I J ) 

111 CONTINUE 

1007 FORMAT (1X»I6) 

WRjTEtLT] ) (NNACC( III »JJ) *JJ = 1 *N$PS) 

IF (MFIN *GE. NSCANS) GO TO 999 
I Y T =MTAR ( 1 ) 

RFADCLT11 ) ( NNACC ( I Y I * JJ ) ♦ JJ = 1 »NSPS) 

MFTN=MFIN+1 
NTFMPsMT A8 ( 1 ) 

IYY=IYYY-1 
DO 121 I = 1 * I YY 
MTAB( I)=MTAB( 1+1 } 

121 CONTINUE 

MTAR ( I YYY ) =NTfMP 
GO TO ?00 
999 CONTINUE 

DO 122 1=2* I YYY 
III =MTAB ( I ) 

DO 112 JJ = 1 .NSPS 

IF (JJ .FQ. NSPS ) GO TO n? 

IF (NNACC( I II » JJ) .LE. NCCNTl GO TO 112 

IF (NNACC (III » JJ+1 ) .LE. NCCNT ) GO TO 112 

IF (NNACCUII.JJ) .LE. 0) GO TO 112 

IF ( NNACC (III* JJ+1 ) .LE. 0) GO TO 112 

IF ( NNACC ( I 1 1 » JJ ) .EQ. NNACC ( I I I » JJ+1 ) ) GO TO 112 

!J*NNACCtIII*JJ) 

Jl=NNACC ( T I I » JJ+1 ) 

IF (JI .GT. 400 .OR. IJ .GT, 400) GO TO 112 
IF (NTRLCJI) .GT. NTBLttJM GO TO 126 
NTRL(IJ)=NTBL(JI ) 

GO TO 11? 

126 CONTINUF 

NTRL( JI)=NTBL( I J) 

112 CONTINUE 

WRlTF(LTl > ( NNACC ( I II *JJ) *JJ = 1 *NSPS) 

122 CONTINUE 
END FILE LTl 
REWIND LTl 
REWIND LT 1 1 
REWIND LT 1 2 

WRlTE(6*l007) (NTRLI I ) *1=1 *400) 

DO 113 1=1*400 

IF ( NTBL ( I ) .FQ. I) GO TO m 
JI = I +1 

IF (NTBL(JI) .NF. I) GO TO 1)4 
NTRL ( JI ) =NTBL ( I ) 

114 CONTINUE 

113 CONTINUE 
11=1 

NTEMPsII 
DO 116 1=2*400 

IF INTbLCII-NTrHI -1 n 117»118*119 
119 IF (NTRL ( I ) .NF. I) GO TO 117 
NTRL ( I “1 ) sNTEMP 


162 



11=11+1 
NTEMPsI I 
60 TO 116 

118 NTBL (1-1 )=NTEMP 
60 TO 116 
117 N=NTRL ( 1 ) 

NTRLU-1 }=NTFMP 
NTFMPsNTRL(N) 

116 CONTINUE 

NTBL(400)=NTEMP 

WR!TF(6*l007) (NTRL( I ) *1=1 *400) 

LWER=1 
LOW=NSTART 
705 CONTINUE 

NUPPER=LWER+120-1 
NHI=L0W+1 ?0-] 

IF (NUPPER *GT. NSPS ) nupper=nsps 

IDIFsNUPPER-LWER+1 

WR I TE ( 6 » 1 006 ) 

1005 FORMAT ( 1 HI ) 

CALL LABEL7(L0W.NHI *1 ) 

00 710 11=1 tNSCANS 

rfad ( LT 1 ) (NNacc( 1 *jj) *JJ=1 * NSPS ) 

DO 115 JJ=1*NSPS 
IB=NNACC( 1 *JJ) 

IF ( IB .LE. 0) GO TO 115 
NNACC(1*JJ)=NTBL{ IB) 

115 CONTINUE 

IF (NFLGXX ,GT, 0) GO TO 127 
WRlTF<LTl2> (NNACCM *J J ) »JJ = 1*NSPS) 

127 CONTINUE 
J!=0 

DO 711 JJ=LWER»NUPPER 
JI=JI+1 

N=NNACC( 1 » JJ)-(NNACCll*JJ)-l >/45*45+2 
I PRT ( J J ) =NSYM ( N ) 

1 PLOT ( J I )=NSYM(N) 

711 CONTINUE 

WRiTf( 6»1005) IT , < 1 PRT ( J J ) . J J=LWpR .NUPP eR 1 
CALL PLTRF7UPL0T, I0IF»NbLK,iNcX* IN C Y,NSTX*NSTY. 
1NCRF*NFLAGX»NFLAG1 ) 

1003 F0RMAT(4X*I6»lH**120Al) 

710 CONTINUE 

REWIND LTl 
NFLAGX=0 
NFLGXX=1 
NFLAGl=0 
. LWER=NUPPER+1 
LOW=NHI +1 

IF (NUPPER .LT, NSPS ) GO TO 705 
570 CONTINUE 

NCCNT=NCNT-1 
END FILE LTl? 

REWIND LT 1 2 

rewind lti 
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IXXX=IXXX-4 
IYYY=IYYY~4 
LT 1 1=13 
RETURN 
END 


SORIGTN BRAVO t SYSUT2*REW 

5JBFTC MOD7B 

SURROUTiNF SfQMRg(NcLUST,KOUNT,!NiTcL ) 

COMMON /LAB1/XRAR<42»12) ,S I GMA ( 42 » 1 2 ) .ROT ( 42 . 1 2 , 1 2 ) 

COMMON /LAR2/X( 12) *aLPhA( 49) ,NsPS »NSC ANS , NCHAN , LT9 » LT 1 0 ♦ LT 1 1 ♦ LT 12 , 
1LT13»LT1 *IXXX»IYYY. 
lNSTART .NSTOP « 

1NBTLG»M0DE»ITYPE*MSFC» I4*NCRE» 

1NSK!P.INCX»INCY,NSTX»NSTY 
DOUBLE PRECISION A ( 1 2 * 1 2 > » EIGEN < 1 2 * 12 ) 

DIMENSION MfRGFC 150) »MPOP( 150) »NFXFC(20) * C ( 42 , 78 ) * B < 12 , 1 2 ) 

dimension com<?4) 

EOUIVALENCE(COM( 1 ) ,NSPS) 

1000 FORMATdX* I6*12F10*?) 

1001 format (1X.4HXPAR ) 

1002 F0RMAT< 1X*16HDID not conveRg ) 

1003 FORMAT < IX » 7HI CLUST= , 16 , 14hMfRgF f ICLUST )= ,16) 

1004 FORMAT (1X,5HRH0= .El 5 . 7 ♦ 5HERR = .E15.7) 

1005 FORMAT (1X.12F10.4) 

1006 FORMAT ( IX ♦ 1 2 1 6 ) 

1007 FORMAT ( IX »23HMeRGI NG WILL TAKE PLACE > 

1 00 A FORMATdH ) 

1009 FORMAT (13H COV. MATRIX ) 

1010 FORMAT ( 1 2H NORM FT GEN ) 

1011 FORMAT (18H P.A. COV. MATRIX ) 

1012 FORMAT < lH *6HASUM= . E 1 5 . 7 ♦ 7HCLUSTER ♦ 1 4 ) 

1013 FORMAT ( IX .28HXRAR ( I » J ) » J=1 » 12 ) » 1=1 »KOUNT ) 

1014 FORMAT < 1X.29HSIGMA ( I ♦ J ) ♦ J = 1 * 12 ) * ! = 1 .KOUNT ) 

1015 FORMAT ( IX *5 5hR0T < I ♦ ICHaN . JchaN ) , JCHAN= 1 . 12) » I CHAN=1 .12 ) ,1 = 1 , KOyNT 
1 ) > 

1016 FORMAT ( 1 X ♦ 1 6 * ( 1 2F1 0*3 ) ) 

1017 F0RMATUH1 ) 
nflg=o 

CZFCH=FLOAT (NCHAN)-?.0 

IF (CZECH .LE. 0.0) CZECH=1.0 

REWIND LT 1 0 

REWIND LT 1 2 

RhO*1 .0/( 10.0**5 ) 

IF ( NSK I P .FQ. 0) GO TO 6 

DO 7 I =1 »NSKI P 

CALL SKRBlN(LTlO,l,NOP) 

7 CONTINUE 

6 CONTINUE 
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DO 5 1CLUST = 1 *NCLUST 
MERGF( ICLUST )=ICLUST 
5 CONTINUE 
IMsNCHAN 
WRJTF(6*1 017) 

IN* IM 
IOPT*l 

DO 10 ICLUST=tNlTCL*NcLUST 
IF (NFLG .GT. 0) GO TO U 
IF ( KOUNT ,GF. 4 2 ) GO TO n 

kount=kount+i 

IFLAGxKOUNT 

CALL FFTCOR ( I FLAG » C* MPOP *NFLG . I N I TCL ) 

WRlTF(6*l 001 ) 

WR I TF < 6 * 1 000 ) I FLAG * ( XRAR ( I FLAG* I ) *1 = 1.12) 

MI * 1 
MJ = 1 2 
MIC* 1 ? 

DO 500 MM = 1 , 1 2 

WRtTf( 6.1000) I FLAG* <C< I FLAG* MR) *MR=Ml ,MJ) 

Ml =MJ+1 
MJsMJ+M!C~MM 
500 CONTINUE 

CALL AMTRX ( IFLAG»XRAR*C*A*NCHAN) 

WRITFI6.1008 ) 

WR I TF ( 6 * 1 009 ) 

WRtTf( 6»1 005 ) ( < A(MT »MJ) »MJ = i ,12 ) .Ml =1,12) 

CALL DJCORI (A»IM,tN,iOPT,RhO,FRR*FIGFN) 
WRlTF(6»lO04) RHO.FRR 

WR]TFt6»l 0051 t ( A ( MI *MJ) *MJ=1,12) .Ml =1*12) 

WR I TF t 6 » 1 008) 

WRlTF(6*l005 ) ( (FIGFNIMI *MJ) *MJ=1 *12) ,Ml =1,12) 

IF (FRR .FQ. 0.0) GO TO 15 
MFRGFUCLUST )=0 
KOUNT=KOUNT-l 
WRITF(6»1002) 

GO TO 10 
15 CONTINUE 

CALL ROTA (I FL AG .ROT *F IGEN *NCHAN ♦ A . SI GMA ) 

MFRGF UCLUST ) =KOUNT 

WRITE ( 6 * 1003) ICLUST*MpRGF( ICLUST ) 

MPOP ( (COUNT ) rMPOP (ICLUST) 

IF (KOUNT •EQ. 1) GO TO 10 

mclust=kount-i 

DO 20 ICHFCK=1 *?0 
NFXFC< ICHFCK)=0 
20 CONTINUE 
MCHFCK=1 

DO 25 JCLUST*) .NCLUST 

IF( MFRGF ( JCLUST ) ,LT .MCHECK )GO TO 25 

MCHFCK=MCHFCK+1 

IF (MFRGFUCLUST) .FO. KOUNT) GO TO 26 
JFLAG=MFRC,F( JCLUST) 

DO 30 ICHAN=l *NCHAN 

X< I CHAN > =XBAR ( JFLAG * I CHAN )-XBAR ( KOUNT ♦ I CHAN ) 
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30 


25 

26 


501 


31 


40 

35 


55 

50 

45 


65 

60 


CONTINUE 

IFLAGaKOUNT 

CALL KCHECK ( I FLAG » ROT » X * SI GMA * ASUM , NCHAN ) 

WR!TE(6*1 008) 

WRlTF(6.1 012>4 SUM»JfLaC, 

IF USUM .GT. CZECH) GO TO 2 5 
IFLAG*JFLAG 

CALL KCHECK I IFLAG » ROT * X . Si GMA ♦ ASUM , NCHAN ) 

WRITE (6*1 000 ) 

WRITE! 6*1 Ol? >ASUM, JFLAG 
IF (ASUM .GT. CZFCH) GO TO 25 
NEXEC(1)=NEXEC(1 )+I 
NSUB=NFXFC« 1 >+l 
NFXEC(NSUB)=JFLAG 


CONTINUE 

IE (NEXEC(l) .EQ. 0) GO TO 10 
DO 501 KK=1 *NSUB 
WRITE(6*1006 ) <K,NfXEC(KK) 
CONTINUE 
MSUB=NEXFC( 1 )+l 
TOTAL=MPOP( KOUNT) 

DO 31 1 RUN = 2 * MSUB 
NSUB=NEXEC( IRUN) 
SUM=MPOP(NSUB) 

TOTAL =TOTAL+SUM 
CONTINUE 


I NUM=0 

DEN=MPOP( KOUNT ) 

DO 35 ICHAN = 1 » NCHAN 

X ( ICHAN ) =XBAR ( KOUNT * I CHAN ) *DFN /TOTAL 
DO 40 JCHAN=ICHAN*NCHAN 
INUM=INUM+1 

B ( ICHAN* JCHAN ) =C ( KOUNT* INUM ) *DEN/TOTAL 

CONTINUE 

CONTINUE 

DO 45 IRUN=?,MSUB 
NSUB=NEXEC( IRUN) 

INUM=0 

DEN=MPOP( NSUB) 

DO 50 I CHAN=1 »NCHAN 

X( ICHAN) =X< ICHAN )+XRAR(NSUP * I CHAN) *DEN/TOTAL 
DO 55 JCHAN=ICHAN»NCHAN 
1 NUM= I NUM+1 

BUCHAN. JCHAN )=P( I CHAN » JCHAN ) +C ( NSUB . I NUM ) *OEN/TOT AL 

CONTINUE 

CONTINUE 

CONTINUE 

DO 60 ICHAN=1 .NCHAN 
DO 65 JCHAN= I CHAN .NCHAN 

A ( ICHAN. JCHAN ) =B( ICHAN. JCHAN )-X( I CHAN ) *X ( JCHAN ) 
AIJCHAN. ICHAN )=A(T CHAN. JCHAN ) 

CONTINUE 
CONTINUE 
WRITF(6»1 00Q) 

WRITE (6. 1005) <<A(M!.MJ) .MJ=1,12) »MI=1.12) 
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CALL DJCOBI (A.IM.IN.IOPT.RhO.ERR.EIGEN) 
WRlTE(6»l004) RHO,ERR 

WRlTF(6»l005) {(A<MT»MJ)tMJ=i»12)*Mi=l*12> 

WRITF(6*T 008 ) 

WRTTF(6*10051 ( <F!GFN(MI.MJ>*MJ=1»12)»MI=1»12) 

IF (FRR .NE. 0.0) GO TO 10 
WRJTF(6*1 007) 

IFLAG*NFXFC(2 ) 

MPOP ( I FLAG ) =TOTAL 
INUMsO 

DO 70 ICHAN=1 »NCHAN 
XRAR( I FLAG* I CHAN ) = X(ICHAN) 

DO 75 JCHAN=1CHAN.NCHAN 
I NUM= INUM+1 

C ( I FLAG * I NUM ) =B ( I CH AN. JCHAN) 

75 CONTINUE 
70 CONTINUE 

call rota ( tflag»Rot»figen*nchan, A*SIGMA) 

DO 80 JCLUST=1 .NCLUST 
DO 85 IRUN=?*MSUB 
NSUB=NEXFC< IRUN) 

IF (MERGFt JCLUST ) ,NE. NSUR) GO TO 85 
MFRGF< JCLUST)=IFLAG 
85 CONTINUE 
80 CONTINUE 

MERGE I ICLUST) = TFLAG 
iFlNFXECm.FO.lJGO TO 94 
ISW=0 
JCHECK=1 

DO 90 JCLUST =1 .NCLUST 
IDUMsMERGE(JCLUST) 

91 IFIMERGE (JCLUST) ,LT.JCHECK)GO TO 90 
IF(MFRGE( JCLUST) .GT.JCHECKlGO TO 92 
IF( ISW.EQ.l )G0 TO 95 

JCHECK= JCHECK +1 

IF( JCHECK. EQ.KOUNT)GO TO 94 

GO TO 90 

92 MFRGF< JCLUST )=MERGF< JCLUST)-1 
ISW=1 

GO TO 91 

93 IF (JCHECK .GT. KOUNT) GO TO 94 
!SW=0 
tNUMaO 

MPOP ( JCHFCK ) sMPOP ( T HUM ) 

DO 95 ICHAN=1 »NCHAN 

XBAR ( JCHECK * ICHAN)=XBAR(IDUM»ICHAN) 

SIGMAIJCHECK* ICHAN)=S1GMA( IDUM. ICHAN) 

DO 100 JCHAN=ICHAN,NCHAN 
1NUM=INUM+1 

C( JCHECK*INUM)=C(IDUM,INUM) 

ROT (JCHECK. I CHAN »JCHAN) =ROT ( I DUM ♦ I CHAN » JCHAN ) 
ROT( JCHFCK, JCHAN. I CHAN )=ROT( IDUM. JCHAN.ICHAN) 
100 CONTI NUF 
95 CONTINUE 

DO 96 LCLUST=1.NCLUST 
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IF(MFRGFfLCLUST) .NF*IDUM)GO TO 96 
MERGE O-CLUST >=JCHECK 
96 CONTINUE 

JCHECK=JCHECK+1 
90 CONTINUE 

94 KOUNT=KOUNT-NFXFC< 1 ) 

10 CONTINUE 

11 CONTINUE 

WRHFU-T9) (COM( I) , j = l ,24) 

WRJTFILT9) { (XraR( I *J) .1 = 1 tlCOUNTi *J»1 ,12) 

wr t tf i LT9 ) ( (Sigma ( i *J) *i=i ,kount ) , j=i , 1 ? > 

WR I T F ( L T 9 ) ( ( (ROT ( I . I CHAN .JCHAN } . 1=1 .KOUNT ) » ICHAN=1 .NChAN ) . 

1 JCHANsl .NCHAN) 

WRITF(6»1 0)3) 

DO 510 1=1 *K0UNT 

WR I Tf ( 6 . 1 000 ) i,(XraR(I»J).J=1»12) 

510 CONTINUE 
WRITE(6*l0l4) 

DO 511 1=1. KOUNT 

WRiTe-(6.lO00) T . (STGMa ( I ♦ J ) . J = 1 . 12 ) 

511 CONTINUE 
WRITF(6»1015) 

DO 512 I =1 » KOUNT 

WRiTF(6»1016) I . { (ROT ( I . I CHAN , JCHAN ) , JCHAN=1 , 12 ) , ICHAN = 1 , 12 ) 

512 CONTINUE 

DO 513 I=1.NCLUST 

IF(MFRGE( I) •GT.KOUNT)GO TO 514 

WRITF(6.615)I .MFRGF( I ) 

515 FORMAT! IX »7HCLUSTfR.I4»1X,5hCLASS. 141 

513 CONTINUE 

514 CONTINUE 

DO 660 1=1 *KOUNT 
DO 620 ICHAN=1 .NCHAN 
DO 610 JCHAN=1*NCHAN 

B( I CHAN. JCHAN) =ROT ( I ♦ JCHAN , I CHAN ) /S I GMA ( I , I CHAN) 

610 CONTINUE 
620 CONTINUE 

DO 650 1CHAN=1 .NCHAN 
DO 640 KCHAN=1 .NCHAN 
SUM*0.0 

DO 630 JCHAN=1 .NCHAN 

SUM=SUM+ROT( I ♦ I CHAN .JCHAN) (JCHAN .KCHAN ) 

630 CONTINUE 

AUCHAN, KCHAN)=SUM 
640 CONTI NUF 
650 CONTINUE 

WR 1 TF ( 6 *600 ) I 

600 FORMAT ( IX . 13HCL ASS FU.lPSE.I4) 

WRITE (6 .1005) ( (AUA.JA) *JA=1*NCHAN)*TA=1 .NCHAN > 

660 CONTINUE 

REWIND LTQ 

RETURN 

END 
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SIBFTF GFTCOR 

SUBROUTINE FETcORt IFLAG»C»NPOP,NFLG*N ) 

COMMON /LAB1/XbaR<42»12) *S IGMa ( 42 * 12 ) *ROT ( 42 * 12 * 12 ) 

COMMON /LAB2/X( 12) » ALPha(^9) *NSPS *N$C ANS * NCHAN *LT9 ♦ LT 10* LT 1 1 *LT12* 
1LT1B*LT1 *1X»IY, 

INSTART *NSTOP* 

1NBTLG»M0DE» ITYPE»MSFC.I4»NCRF* 

1NSKIP*1NCX . INCY»NSTX*NSTY 
DIMENSION NPOP< 1 ) »c(42»78> *NDAT<255 ) 

DATA NCNT/O/ 

1NUM=0 
NFLG1 =0 
NFLG3=0 

DO 5 ICHAN=1 *NCHAN 
XBAR< I FLAG* ICHAN ) = 0.0 
DO 10 jchan=ichan*nchan 
jnum=inum+i 

CUFLAG*INUM) = 0.0 
10 CONTINUE 
5 CONTINUE 
KNT = 0 

AO CONTINUE 

IF (NCNT .GE. NSCANS) GO TO 70 
NFLG?=0 

READtLTl?) (NDAT<JJ),JJ=1*NSPS) 

ncnt =ncnt+i 

NFLGl 

NFLAG?=1 

DO 20 JJ=1*NSPS 

call gft 7 (X(1)»nsps*o»nchan*nscano,ltio*teRP*nflag2* 

INST ART »NbTLG» MODE »NCRE» ITYPF<MSFC ) 

IF (NDAT(JJ) .NE. N) GO TO 30 
KNT=KNT+! 

NFLG2=1 
AI=FLOAT (KNT ) 

I NIJM = 0 

DO 25 ICHAN=1 *NCHAN 

XBAR ( I FLAG* ICHAN > = ( 1 .0-1 .0 /A I )*XBAR ( I FLAG » I CHAN ) +X ( ICHAN ) /A! 

DO 26 JCHAN=ICHAN*NCHAN 
INUM=INUM+1 

C ( I FLAG* I NUM ) = ( 1.0-1,0/AI ) *C ( I FL AG * I NUM ) +X < ICHAN )*X ( JCHAN ) / A I 
26 CONTINUE 
25 CONTINUE 
GO TO 20 
30 CONTINUE 

IF (NFLG2 .FQ, 1 ) GO TO 20 
IF (NDAT(JJ) .NE. N+l ) GOTO 20 
NSAVsNCNT 
NFLG3=1 
20 CONTINUE 

c writf(6»iooo) ncnt,knt*nsav*nbcxup»n*nflg2»nflg3 

1000 FORMAT ( 1 X » 7 T 5 ) 

IF (NFLG? .NF. 0) GO TO 40 
IF (NFLG3 .FO. 0) GO TO 40 

nbckup=ncnt-n$av+i 
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CALL BSRfCD(LTiO,NrcKuP*I4*RF) 
CALL BSRFC0(LT12»NRCKUP»RF ) 

ncnt=nsav-i 

NPOP ( N ) =KNT 

NsN+1 

RFTURN 

70 CONTINUE 

REWIND L T 1 0 
REWIND LT12 
NFLG=1 
RFTURN 
END 


SIRFTF TRXAM 

subroutine AMTRX ( iflag»xbar»c*a»nchan ) 

DIMENSION XRAR<42*1?) *C<42*78) 

DOUBLE PRECISION A ( 1 2 ♦ 1 7 ) 

INUM=0 

DO 1 ICHAN=1 ,NCHAN 
DO 2 JCH A N= I CHAN ♦ NCH AN 
I NUM= I NUM+ 1 

A ( I CHAN ♦ JCH AN ) = C ( I FLAG * I NUM ) -XBAR ( I FLAG* I CHAN ) *XBAR ( I FLAG * JCHAn ) 
At JCHAN* I CHAN > =A (1 CHAN * JCHAN ) 

2 CONTINUE 

1 CONTINUE 

RETURN 
END 


$IBFTC ATOR 

SUBROUTINE ROT At I FL AG » ROT * EIGEN *NCHAN , A * S TGMa ) 
DIMENSION ROT (42*12*1?) 

DIMENSION SIGMa (42*121 
DOUBLE PRECISION AMM2) 

DOUBLF PRECISION FTGEN<12»12) 

DO 1 I CHAN = I * NCH AN 

SIGMA 1 1 FLAG, I CHAN I =A( I CHAN* KHAN) 

DO 2 JCHAN=I *NCHAN 

ROT ( 1 FLAG ♦ I CHAN * JCHAN )=E I GEN (JCHAN* I CHAN) 

2 CONTINUE 

1 CONTINUE 

RETURN 
END 


SIBFTF KCEHCK 

SUBROUTINE KCHFCK ((FLAG* ROT » X * SI GMa » ASUM ,NCHAN ) 

DIMENSION ROT (42*12*12) » SI GMA ( 42 * 1 2 5 *Xf 1) 

ASUMrD.O 

DO 3 ICHAN=1 *NCHAN 
S))MsO.O 

DO 4 JCHAN=1 *NCHAN 

SUM=SUM+ROT ( IFLAG* T CHAN » JCHAN )*X ( JCHAN ) 

4 CONTINUE 
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3 


ASUM=aSUM+SUM*SUM/SIGMa< iflag» ICHAN) 
CONTINUE 

return 

END 


SORIGJN rRaV0.SYSUT2.REW 

IIBFTC M0D7C 

SUBROUTINE CLASS(NCLASS»NTEST,NPASS ) 

COMMON /LAB1/XraRU2»12) » S IGMa ( 42 . 1 2 ) »ROT ( 42 » 1 2 * 1 2 ) 

COMMON /LAR2/X ( 12 ) , ALPHA ( 49 ) , NSPS »N SCANS , NCHAN , LT9 , LTIO »LT1 1 , LT 1 2 » 
1LT13»LT1 ,IX»NDUMMY, 

1 NST ART »NSTOP , 

1 NpTLG i MODE * I TYPF*MSFC , !4 »NcRF* 

1NSKIP,INCX»INCY,NSTX»NSTY 
DIMENSION W ( 1 2 ) » MT AP ( 3 ) 

DIMENSION NDAT(255»3) »PRNT (255) 

DIMENSION COM (?4) 

EOUIVALENCF(COM( 1 ) , N S P S ) 

REWIND LT9 
REWIND LT 1 0 
REWIND LT 1 

Rewind lti2 
rewind lti3 
nelagi=o 
L T 1 = I 

IE (NPASS .Ne. NTFST) LT 1 =LT 1 3 

CZECH=NCHAN 

IE (NSKIP .EO. 0) CO TO 60) 

DO 60? 1=1, NSKIP 

CALL SKRRlN(LTiO,l ,NOP) 

60? CONTINUE 
601 CONTINUE 

READ(LT9) (COMU ) ,T=1,?4) 

Read ( ltq ) ( (xraR ( t »J) .1=1 *nclass) *J= i,i2 j 
RpaD(LT9) C (SIGMA(T*J)*I=1 *NcLaSS) , j=l,l?) 

Read ( LT9 ) ( ( (ROT ( t , i ch an, jch an ) , r = i , nclasS) » ichan=i , nchan » , 

1 Jew AN = 1 , NCHAN ) 

DO 1 1=1,3 
1 MTAB tl ) = T 

DO 10 1END=1 ,NSCANS 

Read ( lt i ? ) (Ndat< i ,1 ) ,1 = 1 ,nspsi 

NELAG2=1 

DO 20 I SUBN = 1 ♦ NSPS 

CALL GET 7 ( X ( 1 ) , N SPS, 0 , NCHAN, NSC ANO, LTIO, IFRR»NfLAG2» 
1NSTART,NRTLG,M0DE»NCRE*ITYPF,MSEC ) 

IE (NDAT ( ISUPN,1 ) .GT. 0 ) NOAT ( I SURN ♦ 1 ) =0 
SMALL=1.75*CZECH 
DO 25 ICLASS=) ,NCLASS 
DO 30 ICHAN=1 ,NCHAN 

W( ICHAN) =X( ICHAN )-XPAR( I CLASS ,TCHAN) 

30 CONTINUE 
A e UM=0, 0 

DO 35 ICHAN=1 , NCHAN 
SUM* 0.0 
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DO 40 JCHAN=1 .NCHAN 

SUM=SUM+Wt JCHAN>*ROT< i class* ICHAN ♦ JCHAN ) 

40 CONTINUE 

41 aSUM=aSUM+SUM*SUM/STGMa ( i class* ICHAN) 

26 CONTTNUr 

IF tASUM .GT. SMALL) GO TO 25 
small=asum 

NdATI ISUpN*1 )=ICLaSS 
25 CONTINUE 
1016 FORMaT(1X*3I6»F10.2) 

20 CONTINUE 

WRlTF(LTl) (NDATt I *1) *1=1 »NSPS) 

10 CONTINUE 

END FILE L T 1 
REWIND LT1 
REWIND LT 1 2 
REWIND LT 1 0 

IE INPASS .Nf. NTFST) GO TO 804 
DO 610 IZ=1*NSC4NS 
I Y=MT AS ( 1 ) 

RFADtLTl ) (NDATt TA.tY)*IA = 1 *NSPS) 
NTFMP=MTAB t 1 ) 

MTAB(1)=MTAR(2) 

MTAPt?)=MTARt2) 

MTAB ( 2 ) =NTEMP 

IE (IZ .LT. 2) GO TO 610 

DO 620 IA=NSTART,N?t0P 

IF ( I A .EQ* 1 ) GO TO 6?0 

IE I I A + 1 .RT. NSTOP ) GO TO 620 

IIY=MTAR< 2 ) 

I MssMTAB ( 1 ) 

I NrMTAB ( 2 ) 

M=NDAT ( I A * I M ) 

N=NDAT( I A — 1 » IN) 

IE t M .NF. N) GO TO 650 
IL=MTAB(2 ) 

L=NDAT(!A»IL ) 

IF (V .NF. L) GO TO 650 
NDATt IA*IIY)=M 
GO TO 620 
660 IMsMTAB(l) 

M-NDATtIA-1 *IM) 

1N=MTAB(5 ) 

N = NDAT( TA-1 * I N ) 

IF (M .NF. N) GO TO 620 
LsNDATM A + i *TM) 

IF CM .NE. L) GO TO 620 
NDATt IA*IIY)=M 
620 CONTINUE 

IF tlZ .LT. 2) GO TO 610 
L=MTARt 1 ) 

WRj Tf t LTl 2 ) tNpATt T *L) . 1=1 .NSPS) 

610 CONT1NUF 

DO 611 1=2*2 
L=MT ARt I ) 
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WRlTF(LTlS) (NoATl !L*L) *TL*l*NSPS) 

611 CONTINUE 
806 CONTINUE 

REWIND LT 1 
REWIND LTP 
REWIND LT18 
REWIND LT1 0 
LOW=NSTART 
LWER=1 

800 CONTINUE 
NHI=L0W+1 ?0-l 
NUPPFR=LWER+1 20-1 

IF (NUPPFR .GT. NSPS ) NUPPFR=NSPS 

IDIFsNUPPER-LWER+1 

WRITF(6»1 007) 

1007 FORMAT ( 1 HI 1 

CALL LABFL7<LOW.NHT *1 ) 

DO 801 11=1 *NSCANS 

READ(LTi-> ) (NDATI JJ,1 ) *JJ = 1 »NSPS) 

DO 808 JJ=LWER. NUPPFR 
IBrNr\AT( JJ*1 ) 

I RND= I B- f IB-1 ) /6 8*68+ 2 
PRNT( JJ)=ALPHA( IRND) 

808 CONTINUE 

CALL PLTBF7(PRNT(LWFR) *IDIF*NBLK» INCX » I NcY *N$TX *NSTY *NcRE ♦ 
1 NFLAGXtNFLAGl ) 

wrjTf(6»i 008) n» (PRnt( JJ) »JJ=lweR*nuppfRi 

1008 FORMAT ( 6X » I 6 ♦ lH* * 1 2 0 A 1 ) 

801 CONTINUE 
REWIND LT 1 8 
NELAG1 =0 
nflagx=o 
lweR=nupper+i 

L0W=NHI+1 

IE (NUPPFR .LT. NSPS ) GO TO 800 

802 CONTINUE 
nscans=nscans-7 
return 

END 
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